Program OVERSPOT
*
c     This routine was originally written in Meteorological Fortran
c     (MeteFor), an extension of Fortran-77, in order to utilize the
c     vector and matrix notation available in MeteFor.  This
c     routine is also maintained in MeteFor.  Some of the original
c     MeteFor code may appear in statements which have been commented
c     out.  The original MeteFor source (suffix .hlf) is more
c     readable and self-documenting.  See the document MeteFor.doc.
*
*     Inputs:
*
c     Satellite: (AQUA, TERRA, NOAA18, etc)
c     Starting date: yymmdd
c     Ending data:   yymmdd
c     Geodetic Latitude of Ground Point
c     Longitude of ground point
c     Target Diameter in kilometers
*
	implicit double precision(d)
	parameter(d0=0.d+0, dtintv=48.d+0/1440.d+0)
	character*4 cmode
	character*8 csat
*     vector (f)vcoord,(f)vnadir,(f)vec4,(f)vnorax
	real vq41(3),vq42(3)
	real  vspot(3),vwb(3),vsat(3),vnad(3),vdl(3)
cxr   numarg,fltarg,dtmjdn,vec4,vcoord,vnadir,vnorax,invnor,angbtw
cxr   tcivil
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	call getarg(1,csat)
	jd1 = argnum(2)
	jd2 = argnum(3)
	gdlat = fltarg(4)
	glong = fltarg(5)
	sepmax = fltarg(6)
	dtgo = dtmjdn(jd1/10000, mod(jd1/100,100), mod(jd1,100),0,0,0)
	dtstop = dtmjdn(jd2/10000, mod(jd2/100,100), mod(jd2,100),0,0,0)
	open(30, file='OVERSPOT', form='formatted', status='unknown',
     x  access='sequential')
	write(30, '(1h ,a8)') csat
	re = erad84(gdlat)
*
*     vspot = vcoord(d0, vec4(gdlat,glong,re), 'DLT')
	call vec4(gdlat,glong,re,vq41)
	call vcoord(d0,vq41,'DLT',vq42)
	vspot(1) = vq42(1)
	vspot(2) = vq42(2)
	vspot(3) = vq42(3)
	nc = 0
*
	do   2 dt=dtgo,dtstop,dtintv
*     nc = & + 1
	nc=nc+1
*
	if(mod(nc,100).eq.0)then
*     vwb=vnorax(d0,'        ',wb,wb,wb)
	call vnorax(d0,'        ',wb,wb,wb,vwb)
	end if
*
	call invnor(dt,csat,  gdlat,glong,
     x  dtpass,slant,elev,scan,sclat,sclong,cdt,cmode)
	if(dtpass.lt.dt)go to 2
*
*     vsat = vcoord(d0, vec4(sclat,sclong,cdt), 'LLT')
	call vec4(sclat,sclong,cdt,vq41)
	call vcoord(d0,vq41,'LLT',vq42)
	vsat(1) = vq42(1)
	vsat(2) = vq42(2)
	vsat(3) = vq42(3)
*
*     vnad = vnadir(vsat)
	call vnadir(vsat,vnad)
*
*     sep = 111.12*angbtw(vnad,vspot)
	rq41=angbtw(vnad,vspot)
	sep=111.12*rq41
	if(sep.gt.sepmax)go to 2
*
	call tcivil(dtpass, jy,jm,jd,kh,km,sec,nth)
*     anad = angbtw(vnad-vsat, vspot-vsat)
	vq41(1) = vnad(1)-vsat(1)
	vq41(2) = vnad(2)-vsat(2)
	vq41(3) = vnad(3)-vsat(3)
	vq42(1) = vspot(1)-vsat(1)
	vq42(2) = vspot(2)-vsat(2)
	vq42(3) = vspot(3)-vsat(3)
	rq41=angbtw(vq41,vq42)
	anad=rq41
*
*     vdl = vcoord(d0, vnad, 'TDL')
	call vcoord(d0,vnad,'TDL',vdl)
*
	write(30, '(1h , i5,4(1x,i2.2),f6.1, 3f9.3,f9.2)') jy,jm,jd,kh,km,sec,
     x  vdl(1),vdl(2),sep,anad
	write(*, '(1h , i5,4(1x,i2.2),f6.1, 3f9.3,f9.2)') jy,jm,jd,kh,km,sec,
     x  vdl(1),vdl(2),sep,anad
    2 continue
*
	close(30)
	stop
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function ARGNUM(jpos)
*
*     15 Jun 06 - uses new nscan
	character*24 cin
cxr   nscan
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	call getarg(jpos, cin)
	kol = 1
	num = nscan(kol,cin)
*     num may be negative
	argnum = num
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function FLTARG(jpos)
*
c     To fetch the jpos-th input argument (like IPP in McIdas)
c     as a REAL
	Character*1 cinput(80), ch
	Character*80 cin
	equivalence(cin,cinput(1))
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	cin  = ' '
	call getarg(jpos, cin)
	kol = 1
	jsign = 1
	nfrac = 0
	kpt = 0
	mag = 0
*
	if(cinput(1) == '-') then
	kol = 2
	jsign = -1
	end if
*
	do n = 1,32
	 ch = cinput(kol)
	 if(ch == ' ') exit
	 if(ch == '.') kpt = 1
	 ich = ichar(ch)
*
	 if(ich>=48 .and. ich<=57) then
	  mag = 10*mag + ich - 48
	  nfrac = nfrac + kpt
	 end if
*
	 kol = kol + 1
	end do
*
	flt = jsign*mag
	if(nfrac > 0) flt = flt / (10. ** nfrac)
	fltarg = flt
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Double Precision Function DTMJDN(myrx,mon,mday,jhour,jmin,jsec)
*
c     To compute the Modified Julian Day Number (MJDN), based on
c     the epoch of 00z 1 January 1970:
*
*     The Modified Julian Day Number is the Julian Day Number
c     of a given date, with the JDN of 12z 1 Jan 1970 subtracted.
*
*      The last 3 arguments may be either REAL*4 or INTEGER*4.
*
*     This routine should not be used for dates earlier than
c     1 January 1970.
*
	Implicit None
	Double Precision djul,dhour,dmin,dsec
	Real hour,fmin,sec
	Integer khour,kmin,ksec,imax,myr,myrx,jhour,jmin,jsec,kalday
	Integer mday,jday,month,mon,jcode,m100,mod,m400,m4,m,ndm
	Integer julday
*
	equivalence(khour,hour),(kmin,fmin),(ksec,sec)
	save
	data imax/16777216/
cxr   julday
c     xxxxxxxxxxxxxxxxxxxxxxxxxxx
	myr = myrx
	khour = jhour
	kmin = jmin
	ksec = jsec
	if(khour .lt. imax)hour = khour
	if(kmin .lt. imax)fmin = kmin
	if(ksec .lt. imax)sec = ksec
	kalday = mday
	jday = kalday
	month = mon
	if(mon .ne. 0) go to 10
*
	jcode = 1115212
	m100 = mod(myr,100)
	m400 = mod(myr,400)
	m4 = mod(myr,4)
	if(m400.eq.0 .or. (m4.eq.0 .and. m100.ne.0))jcode = 1115208
*
	do 12 m = 1,12
	month = m
	ndm = 31 - mod(jcode,4)
	jday = jday - ndm
	if(jday.le.0) go to 20
*
12    jcode = jcode/4
*
20    kalday = jday + ndm
*
	if(month.gt.12.or.kalday.gt.31)then
	write(*, '('' DTMJDN: month & day:'', 2i8)') month,kalday
	stop 99
	end if
*
10    djul = julday(myr, month, kalday) - 2440588
	dhour = hour
	dmin = fmin
	dsec = sec
	djul = djul + (dhour + dmin/60.d+0 + dsec/3600.d+0)/24.d+0
	dtmjdn = djul - .5d+0
	if(djul .le. 0.d+0)dtmjdn = 0.d+0
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
c     Vector*4 Function VEC4(x,y,z)
	Subroutine VEC4(x,y,z,  vec)
c     To return a vector whose components are the three given
c     values x,y,z
*
	real vec(3)
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	vec(1) = x
	vec(2) = y
	vec(3) = z
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Subroutine VCOORD(dt,vecin,cmode,  vecout)
*
cxr   sine,cosine,geodet,geocen,arctan,verneq
! -----------------------------------------------------------------
! DOCUMENTATION:
*
! Cooperative Institute for Meteorological Satellite Studies
! University of Wisconsin
! 1225 W Dayton Street
! Madison, Wisconsin  53706      USA
*
! Name: VCOORD; a subroutine-type subprogram
!     This is a vector-valued function under MeteFor.
*
! Inputs: Time, real(double), in Julian Day Number
!         or Modified Julian Day Number (MJDN)
!     a real(single) 3-dimensional vector in some coordinate system
!     a character*3 variable stating the desired transformation
!     12 Jan 00: Modified Julian Day Number (from epoch of
!     1 Jan 70) is also accepted.
*
! Output: a 3-dimensional real(single) transformed vector
*
! Usage:
!     To inter-convert among various coordinate systems:
!      cmode = LLC         Geocentric Lat/Long to Celestial
!      cmode = LLT         Geocentric Lat/Long to Terrestrial
!      cmode = CT          Celestial to Terrestrial
!      cmode = DLC         Geodetic Lat/Long to Celestial
!      cmode = DLT         Geodetic Lat/Long to Terrestrial
*
!     and the inverses of these operations:  (CLL, TLL, TC, DLT, TDL)
*
!     Note regarding time:  An input value for time (dt) is required
!     only if the celestial coordinate system is involved, either as
!     input or output coordinates.  If the transformation is
!     strictly between Lat/Long/CD on the one hand, and terrestrial
!     on the other, a dummy value of time (e.g. zero) can be used.
*
! Modules:  verneq, sine, cosine, arctan
*
! Comment:  The three coordinate systems treated by this routine are
!   defined as follows:
*
!   Celestial: The x-coordinate is directed from the center of
!   the earth toward the vernal equinox.  The y-coordinate is
!   directed 90 deg east of the x-coordinate, and the z-coordinate
!   is directed from the earth's center toward the north pole.
!   The three coordinates form a right-hand orthonormal system.
*
!   Terrestrial: The x-coordinate is directed from the earth's center
!   toward the Greenwich meridian; the y-coordinate is directed
!   90 deg east of the x-coordinate; the z coordinate is directed
!   from the earth's center toward the north pole.  The terrestrial
!   coordinates are a right-hand orthonormal system.
*
!   Latitude/Longitude:  The three components are latitude in
!   degrees, positive north; longitude in degrees, positive east;
!   and radius vector in arbitrary units from the earth's center
!   to a given point.  These coordinates are not orthogonal
!   vectors.  The coordinates are vectors only in the sense that
!   they are three ordered real numbers.  The latitude may be either
!   geocentric or geodetic, depending on the conversion mode.
*
! Programmer: F W Nagle
!  ---------------------------------------------------------------------
	double precision dt,deqr2,deqrad
	character*3 cmodes(10),cmode
	dimension vecin(3),vecout(3),vecinp(3)
	save
	data cmodes/'LLC', 'LLT', 'CT ', 'CLL', 'TLL', 'TC ',
     x  'DLC', 'DLT', 'CDL', 'TDL'/
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	kdetic = 0
	lm = 0
*
	do   2 j=1,10
	lm = j
	if(cmodes(j)(1:3).eq.cmode(1:3)) go to 8
    2 continue
*
	write(*, '('' VCOORD: Unrecognized Transformation: '',
     1  a3)') CMODE
	return
8     jtest = (lm-1)*(lm-3)*(lm-4)*(lm-6)*(lm-7)*(lm-9)
	velong= 0.
*
	if(jtest.eq.0)then
	if(dt.gt.2.6d+6.or.dt.lt.0.d+0) go to 920
	if(dt.gt.100000.d+0.and.dt.lt.2.4d+6) go to 920
	  velong= verneq(dt)
	end if
*
*
	do   4 j=1,3
	vecinp(j) = vecin(j)
    4 continue
*
	if(lm.eq.7.or.lm.eq.8)then
	lm=lm-6
	vecinp(1)=geocen(vecinp(1))
	end if
*
	if(lm.ge.9)then
*     lm = & - 5
	lm=lm-5
	kdetic = 1
	end if
*
	if(lm.ge.4)lm=3-lm
	if(iabs(lm).gt.1) go to 50
	if(lm.lt.0) go to 20
*
10    if(abs(vecinp(1)).gt.90.)then
	write(*, '('' VCOORD: Bad input latitude: '', f12.2)')
     1  vecinp(1)
	go to 922
	end if
*
	deqrad = vecinp(3) * cosine(vecinp(1))
	clong = vecinp(2) - velong
	vecout(1) = deqrad * cosine(clong)
	vecout(2) = deqrad * sine(clong)
	vecout(3) = vecinp(3) * sine(vecinp(1))
	return
*
20    clong = arctan(vecinp(1), vecinp(2))
	deqr2 = dble(vecinp(1))**2 + dble(vecinp(2))**2
	eqr = dsqrt(deqr2)
	vecout(1) = arctan(eqr, vecinp(3))
	vecout(2) = clong + velong
	if(vecout(2).lt.-180.)vecout(2) = vecout(2)+360.
	if(vecout(2).gt.180.)vecout(2) = vecout(2)-360.
	vecout(3) = dsqrt(deqr2 + dble(vecinp(3))**2)
	if(kdetic.gt.0)vecout(1)=geodet(vecout(1))
	return
*
50    if(iabs(lm).gt.2) go to 100
	if(lm.lt.0) go to 20
	go to 10
*
100   if(lm.lt.0) go to 110
*
102   vecout(3) = vecinp(3)
	clong = arctan(vecinp(1), vecinp(2))
	glong = clong + velong
	deqrad = dsqrt(dble(vecinp(1))**2 + dble(vecinp(2))**2)
	vecout(1) = deqrad * cosine(glong)
	vecout(2) = deqrad * sine(glong)
	return
*
110   velong = -velong
	go to 102
*
920   write(*, '('' VCOORD: The input time is wrong: '', 1pd20.10)')
     1  dt
*
922   vecout(1) = 0.
	vecout(2) = 0.
	vecout(3) = 0.
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Subroutine VNADIR(vsat,vqadir)
*
cxr   vcoord,vec4,arctan,geodet,erad84
*     Given a celestial or terrestrial satellite position vector,
c     this routine computes in the same coordinate system the
c     position vector of the point on the ground at the satellite
c     nadir, taking into account the ellipsoidal shape of the
c     earth.   This point in general is not at the geocentric
c     latitude of the satellite.
*
	implicit double precision(d)
	parameter (smajx=6378.388, sminx=6356.911946,
     1  b=sminx/smajx, b2=b*b, d0=0.d+0, epsi=.000001)
	integer jqmax/z'7fffffff'/
	real vqadir(3)
*     vector (f)vcoord,(f)vec4
	real vq41(3),vq42(3)
	real  vbl(3),vsat(3),vll(3),vs(3),vtan(3),vg(3),view(3)
	save
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
*     vbl = vsat
	vbl(1) = vsat(1)
	vbl(2) = vsat(2)
	vbl(3) = vsat(3)
*
*     vll = vcoord(d0, vbl, 'TLL')
	call vcoord(d0,vbl,'TLL',vq41)
	vll(1) = vq41(1)
	vll(2) = vq41(2)
	vll(3) = vq41(3)
*
	gclat = vll(1)
	kase = 0
	if(abs(gclat).lt.45.)kase=1
	if(kase.eq.0)then
*
*      vs = vbl/smajx
	vq41(1) = vbl(1)/smajx
	vq41(2) = vbl(2)/smajx
	vq41(3) = vbl(3)/smajx
	vs(1) = vq41(1)
	vs(2) = vq41(2)
	vs(3) = vq41(3)
*
*      vs = vec4(sqrt(vs(1)**2 + vs(2)**2), abs(vs(3)), 0.)
	rq41=vs(1)
	rq42=vs(2)
	rq43=rq41**2
	rq44=rq42**2
	rq45=rq43+rq44
	call vec4(sqrt(rq45),abs(vs(3)),0.,vq41)
	vs(1) = vq41(1)
	vs(2) = vq41(2)
	vs(3) = vq41(3)
*
	 c = b
	 c2 = b2
	   else
*
*      vs = vbl/sminx
	vq41(1) = vbl(1)/sminx
	vq41(2) = vbl(2)/sminx
	vq41(3) = vbl(3)/sminx
	vs(1) = vq41(1)
	vs(2) = vq41(2)
	vs(3) = vq41(3)
*
*      vs = vec4( abs(vs(3)), sqrt(vs(1)**2 + vs(2)**2), 0.)
	rq41=vs(1)
	rq42=vs(2)
	rq43=rq41**2
	rq44=rq42**2
	rq45=rq43+rq44
	call vec4(abs(vs(3)),sqrt(rq45),0.,vq41)
	vs(1) = vq41(1)
	vs(2) = vq41(2)
	vs(3) = vq41(3)
*
	 c = 1./b
	 c2 = 1./b2
	end if
	x = vs(1)/sqrt(vs(1)*vs(1) + vs(2)*vs(2))
	xl = x + .01
	sdotj = vs(2)
	vtan(1) = 1.
	vtan(3) = 0.
*
*     do until(abs(x-xl) < epsi)
	do   2  kq0001 = 1,jqmax
	if(abs(x-xl).lt.epsi)go to 4
*
	 x2 = x*x
	 radic = sqrt(1. - x2)
*
*      vg = vec4(x, c*radic, 0.)
	rq41=c*radic
	call vec4(x,rq41,0.,vq41)
	vg(1) = vq41(1)
	vg(2) = vq41(2)
	vg(3) = vq41(3)
*
*
*      view = vs - vg
	vq41(1) = vs(1)-vg(1)
	vq41(2) = vs(2)-vg(2)
	vq41(3) = vs(3)-vg(3)
	view(1) = vq41(1)
	view(2) = vq41(2)
	view(3) = vq41(3)
*
	 vtan(2) = -x*c/radic
*
*      prod =  view*vtan
	rq41=view(1)*vtan(1)+view(2)*vtan(2)+view(3)*vtan(3)
	prod=rq41
*
*      if(abs(prod) <= 0.) exit
	if(abs(prod).le.0.)then
	go to     4
	end if
	 pp = -c*(x2/(radic**3) + radic)*sdotj - 1. + c2
	 xl = x
	 x = x - prod/pp
    2 continue
*
4     y = c*sqrt(1. - x*x)
	if(kase.eq.0)then
	 elat = sign(1.,gclat)*arctan(x,y)
	   else
	 elat = sign(1.,gclat)*arctan(y,x)
	end if
	rad = erad84(geodet(elat))
*
*     vnadir = vcoord(d0, vec4(elat, vll(2), rad), 'LLT')
	call vec4(elat,vll(2),rad,vq41)
	call vcoord(d0,vq41,'LLT',vq42)
	vqadir(1) = vq42(1)
	vqadir(2) = vq42(2)
	vqadir(3) = vq42(3)
*
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
*     VECTOR FUNCTION VNORAX(DT,CSATX,  SLAT,SLONG,CDT)
	SUBROUTINE VNORAX(DT,CSATX,SLAT,SLONG,CDT,VQORAX)
*
c     This routine was originally written in Meteorological Fortran
c     (MeteFor), an extension of Fortran-77, in order to utilize the
c     vector and matrix notation available in MeteFor.  This
c     routine is also maintained in MeteFor.  Some of the original
c     MeteFor code may appear in statements which have been commented
c     out.  The original MeteFor source (suffix .hlf) is more
c     readable and self-documenting.  See the document MeteFor.doc.
*
c     ------------------------------------------------------------
*    Orbital parameters (2-line elements) are available
*    through the web site:  celestrak.com/pub/elements/weather.txt
*
*     1 Sep 07; in util.a
c     21 Jun 07; The common block /tle/ is now used to contain csats and
c      parms, to allow internal parameter updates with INSTOP.
c     24 Oct 06; allows 24 vice 16 satellites, used in NORWATCH
*     17 May 01; initialization of common blocks transferred to
c       INPARM, q.v.
*     16 ===> 8 as maximum num of satellites
*
!     FUNCTION:
*
!     This routine computes the past or future position of
!     any of several satellites at a given time, returning the
!     results as latitude/longitude, or as a 3-dimensional
!     vector in the celestial coordinate system.   The calculation
!     is based on the NORAD-developed routine SGP8.
!     A detailed document in Acrobat format exists describing the
!     NORAD-developed routines (spacetrk.pdf).
*
!     INPUTS:
*
!     dt (Double Precision or Real(double)):
!     The time for which a predicted satellite position is
!     wanted, expressed in Modified Julian Day Number (MJDN).
!     For an explanation of MJDN, see the ASCII document
!     Julian.Day.Number, and the routine dtmjdn.f
*
!     csatx (character*8):
!     The identifier for the satellite (e.g. NOAA11, DMSP14, etc)
*
!     OUTPUTS:
*
!    slat(real or real(single)):
!    satellite geocentric latitude
*
!    slong (real or real(single)):
!    satellite longitude
*
!    cdt (real or real(single)):
!    central distance, measured from the center of the earth
!    to the satellite, in kilometers.
*
!    VQORAX (Vector*8 or real(single), dimension=3):
!    The vector position of the satellite in celestial coordinates.
*
!     Note: The needed orbital parameters must exist in the NORAD
!     two-line element format in the file TLINSET.
*
!     PROGRAMMER: F W Nagle, SSEC
c     ---------------------------------------------------
	implicit double precision(d)
	parameter (nmax=24)
*
	character*8 csat,csatx,cslast,csats
	character*132 caplet
	real de2ra,buf(13)
	real vqorax(3)
*     vector (f)vcoord
	real vq41(3)
	real  vsat(3),vll(3)
	dimension csats(nmax),parms(13,nmax)
*
	common/e1/dtelem,ds50,xmo,xnodeo,omegao,eo,xincl,xno,xndt2o,xndd60
     x,bstar,x,y,z,xdot,ydot,zdot
	common/tle/csats,parms
	common/c1/ck2,ck4,e6a,qoms2t,s,tothrd,xj3,xke,eqrad,xmnpda,ae
	common/c2/de2ra,pi,pio2,twopi,x3pio2
	equivalence(buf(1),dtelem)
	save
*
	data init/0/
	data csats/ 24*' ' /, cslast/' '/
	data xj2,xj4,qo,so/1.082616e-3, -1.65597e-6,120., 78./
cxr   caplet,inparm,sgp8,vcoord
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if(init.eq.0)instal=0
	jsat = 0
	csat = csatx
*     csat = caplet(&,8)
	csat=caplet(csat,8)
	iflag = 0
*
	if(csat.eq.' ')then
*
	do   2 j=1,nmax
2     csats(j) = ' '
*
*      vnorax = 0.
	vqorax(1)=0.
	vqorax(2)=0.
	vqorax(3)=0.
	 cslast = ' '
	 init = 0
	 instal = 0
	 go to 102
*
	end if
*
	if(csat.ne.cslast)then
	 iflag = 1
	 cslast = csat
	end if
*
	do   4 n=1,nmax
*
	if(csats(n).eq.csat)then
	  jsat = n
	  kfetch = 0
	go to     6
	end if
*
    4 continue
    6 continue
*
	if(jsat.eq.0)then
	 kfetch = 1
*
	do   8 i=1,nmax
	if(csats(i).ne.' ')go to 8
*
	  jsat = i
	  csats(i) = csat
	go to    10
    8 continue
*
   10 continue
*
	if(jsat.eq.0)then
	write(*, '('' VNORAX: There is no slot available for this '',
     1  ''satellite: '', a8)') CSAT
*
*     vnorax = 0.
	vqorax(1)=0.
	vqorax(2)=0.
	vqorax(3)=0.
	go to 102
	end if
*
	end if
*
	call tcivil(dt, ny,nm,nd,km,km,sec,nth)
*
	if(kfetch.gt.0)then
	 jymd = 10000*ny + 100*nm + nd
*
	if(inparm(jymd,csat).ne.0)then
	 write(*, '('' VNORAX: Error return from INPARM '')')
*
*      vnorax = 0.
	vqorax(1)=0.
	vqorax(2)=0.
	vqorax(3)=0.
	 go to 102
*
	end if
*
	do  12 i=1,13
12    parms(i,jsat) = buf(i)
*
	end if
*
	do  14 i=1,13
14    buf(i) = parms(i,jsat)
*
	 iflag = 1
	 instal = jsat
	 jte = dtelem
	 dtf = dtelem - dflt(jte)
	 dtepoc = dtmjdn(jte/1000, 0, mod(jte,1000), 0,0,0) + dtf
	tsince = 1440.d+0 * (dt - dtepoc)
	call sgp8(iflag,tsince)
	xdot = xdot * eqrad / 60.
	ydot = ydot * eqrad / 60.
	zdot = zdot * eqrad / 60.
	fac = eqrad/ae
	vsat(1) = fac*x
	vsat(2) = fac*y
	vsat(3) = fac*z
*
*     vll = vcoord(dt, vsat, 'CLL')
	call vcoord(dt,vsat,'CLL',vll)
	slat = vll(1)
	slong = vll(2)
	cdt = vll(3)
*
*     vnorax = vsat
	vqorax(1) = vsat(1)
	vqorax(2) = vsat(2)
	vqorax(3) = vsat(3)
	init = 1
*
102   return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	SUBROUTINE INVNOR(DT0,CSAT, GDLAT,GLONG,
     1  dtpass,slant,elev,scan,sclat,sclong,cdt,cmode)
c
c     This routine was originally written in Meteorological Fortran
c     (MeteFor), an extension of Fortran-77, in order to utilize the
c     vector and matrix notation available in MeteFor.  This
c     routine is also maintained in MeteFor.  Some of the original
c     MeteFor code may appear in statements which have been commented
c     out.  The original MeteFor source (suffix .hlf) is more
c     readable and self-documenting.  See the document MeteFor.doc.
c
cxr   vnorax,dtmjdn,vnadir,vcoord,rvmag4,angbtw,unit4,vconic,tcivil
c     --------------------------------------------------------------
*     Generic Inverse Navigation:  This routine computes the time
c     of transit of a satellite with respect to a given ground point.
c     "Time of Transit" means the moment at which a satellite views
c     a given ground point directly to its right or left.
*
c     Inputs:
c     dt0 is the Modified Julian Day Number (MJDN) at some time
c       within half an orbit, forward or backward, of the unknown
*       time of transit.
c       The epoch of MJDN is 00Z 1 Jan 1970.
*     Longitude (grlong) is positive EAST of Greenwich
*     Latitude (gdlat) is geodetic.
*
*     Outputs:
*     dtpass - time of transit (Double Precision) in MJDN.
*     slant range in kilometers at transit
*     elevation angle in deg as seen from the ground
*     nadir angle of ground point as seen from satellite, negative
c       being to the left of track
*     satellite geocentric latitude in deg at transit
*     satellite longitude
*     central distance (cd) in km from earth center to satellite
*     mode (character*4), Asc or Desc (ascending or descending path)
c     ----------------------------------------------------------------
	implicit double precision(d)
	parameter(d0=0.d+0, dt1min=1.d+0/1440.d+0)
	character*4 cmode
	character*8 csat
	integer jqmax/z'7fffffff'/
*     vector (f)vnorax,(f)vcoord,(f)unit4,(f)vec4,(f)vconic,
*    1  (f)vnadir
*     vector uk/0.,0.,1./
	real uk(3)/0.,0.,1./
	real vq41(3),vq42(3),vq43(3)
	real  vtgr(3),vsat(3),uop(3),ugc(3),uxc(3),vcgr(3),vz(3),view(3),u
     xopr(3),vslant(3),vll(3)
	dimension xdum(16),veloc(3)
	common/e1/xdum,veloc
	equivalence(ecc,xdum(8)),(xmm,xdum(10))
	equivalence(xincl,xdum(9))
	equivalence(dtelem,xdum(1))
	save
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	gclat = geocen(gdlat)
	re = earrad(gclat)
*
*     vtgr = vcoord(d0, vec4(gclat,glong,re), 'LLT')
	call vec4(gclat,glong,re,vq41)
	call vcoord(d0,vq41,'LLT',vq42)
	vtgr(1) = vq42(1)
	vtgr(2) = vq42(2)
	vtgr(3) = vq42(3)
*
*     vsat = vnorax(dt0,csat, slat,slong,cd)
	call vnorax(dt0,csat,slat,slong,cd,vq41)
	vsat(1) = vq41(1)
	vsat(2) = vq41(2)
	vsat(3) = vq41(3)
*
*     if(rvmag4(vsat) <= 0.) then
	rq41=rvmag4(vsat)
*
	if(rq41.le.0.)then
	cdt = 0.
	write(*, '('' INVNOR: VSAT has zero magnitude'')')
	stop 99
	end if
*
	dmm = 1440.*xmm/6.2831853
	dincl = 57.2957795*xincl
	deccen = ecc
	sma = 1.852*smaxis(dmm, deccen, dincl)
	permin = (sma/330.9982746) ** 1.5
	jdelem = dtelem
	dtepoc = dtmjdn(jdelem/1000, 0, mod(jdelem,1000),0,0,0)
     1  +  dtelem - dflt(jdelem)
	dtepsi = .1d+0 + 2.3d+0*dabs(dt0-dtepoc)/1000.d+0
*     dtepsi = &/86400.d+0
	dtepsi=dtepsi/86400.d+0
*
*     uop = unit4(vsat <x> veloc)
	vq41(1) = vsat(2)*veloc(3) - veloc(2)*vsat(3)
	vq41(2) = vsat(3)*veloc(1) - veloc(3)*vsat(1)
	vq41(3) = vsat(1)*veloc(2) - veloc(1)*vsat(2)
	call unit4(vq41,vq42)
	uop(1) = vq42(1)
	uop(2) = vq42(2)
	uop(3) = vq42(3)
*
*
*     ugc = vcoord(dt0, vtgr, 'TC ')/re
	call vcoord(dt0,vtgr,'TC ',vq41)
	vq42(1) = vq41(1)/re
	vq42(2) = vq41(2)/re
	vq42(3) = vq41(3)/re
	ugc(1) = vq42(1)
	ugc(2) = vq42(2)
	ugc(3) = vq42(3)
*
*     uxc = unit4(uop <x> (ugc <x> uop))
	vq41(1) = ugc(2)*uop(3) - uop(2)*ugc(3)
	vq41(2) = ugc(3)*uop(1) - uop(3)*ugc(1)
	vq41(3) = ugc(1)*uop(2) - uop(1)*ugc(2)
	vq42(1) = uop(2)*vq41(3) - vq41(2)*uop(3)
	vq42(2) = uop(3)*vq41(1) - vq41(3)*uop(1)
	vq42(3) = uop(1)*vq41(2) - vq41(1)*uop(2)
	call unit4(vq42,vq43)
	uxc(1) = vq43(1)
	uxc(2) = vq43(2)
	uxc(3) = vq43(3)
*
*     ang = sign( angbtw(vsat,uxc), vsat<x>uxc*uop)
	vq41(1) = vsat(2)*uxc(3) - uxc(2)*vsat(3)
	vq41(2) = vsat(3)*uxc(1) - uxc(3)*vsat(1)
	vq41(3) = vsat(1)*uxc(2) - uxc(1)*vsat(2)
	rq41=vq41(1)*uop(1)+vq41(2)*uop(2)+vq41(3)*uop(3)
	rq42=sign(angbtw(vsat,uxc),rq41)
	ang=rq42
	gammax = arccos( re/cd)
*
*     angb= angbtw(ugc,uxc)
	rq41=angbtw(ugc,uxc)
	angb=rq41
*
	if(angb.ge.gammax)then
	dtpass = -1.d+0
	go to 102
	end if
*
	togo = ang*permin/360.
	dtl = dt0 + togo/1440.
*
*     vcgr = vcoord(dtl, vtgr, 'TC  ')
	call vcoord(dtl,vtgr,'TC  ',vq41)
	vcgr(1) = vq41(1)
	vcgr(2) = vq41(2)
	vcgr(3) = vq41(3)
*
*     vsat = vnorax(dtl,csat, slat,slong,cd)
	call vnorax(dtl,csat,slat,slong,cd,vq41)
	vsat(1) = vq41(1)
	vsat(2) = vq41(2)
	vsat(3) = vq41(3)
*
*     vz = vsat - vnadir(vsat)
	call vnadir(vsat,vq41)
	vq42(1) = vsat(1)-vq41(1)
	vq42(2) = vsat(2)-vq41(2)
	vq42(3) = vsat(3)-vq41(3)
	vz(1) = vq42(1)
	vz(2) = vq42(2)
	vz(3) = vq42(3)
*
*     view = vcgr - vsat
	vq41(1) = vcgr(1)-vsat(1)
	vq41(2) = vcgr(2)-vsat(2)
	vq41(3) = vcgr(3)-vsat(3)
	view(1) = vq41(1)
	view(2) = vq41(2)
	view(3) = vq41(3)
*
*     spl = view * (uop <x> vz)
	vq41(1) = uop(2)*vz(3) - vz(2)*uop(3)
	vq41(2) = uop(3)*vz(1) - vz(3)*uop(1)
	vq41(3) = uop(1)*vz(2) - vz(1)*uop(2)
	rq41=view(1)*vq41(1)+view(2)*vq41(2)+view(3)*vq41(3)
	spl=rq41
	dt = dtl + dt1min
	nits = 0
c
*     do until (dabs(dt-dtl) < dtepsi)
	do   2  kq0001 = 1,jqmax
*
	if(dabs(dt-dtl).lt.dtepsi)go to 4
*
*      vsat = vnorax(dt,csat, slat,slong,cd)
	call vnorax(dt,csat,slat,slong,cd,vq41)
	vsat(1) = vq41(1)
	vsat(2) = vq41(2)
	vsat(3) = vq41(3)
*
*      vcgr = vcoord(dt, vtgr, 'TC  ')
	call vcoord(dt,vtgr,'TC  ',vq41)
	vcgr(1) = vq41(1)
	vcgr(2) = vq41(2)
	vcgr(3) = vq41(3)
*
*      uop = unit4(vsat<x>veloc)
	vq41(1) = vsat(2)*veloc(3) - veloc(2)*vsat(3)
	vq41(2) = vsat(3)*veloc(1) - veloc(3)*vsat(1)
	vq41(3) = vsat(1)*veloc(2) - veloc(1)*vsat(2)
	call unit4(vq41,vq42)
	uop(1) = vq42(1)
	uop(2) = vq42(2)
	uop(3) = vq42(3)
*
*      uopr = uop
	uopr(1) = uop(1)
	uopr(2) = uop(2)
	uopr(3) = uop(3)
*
*      vz = vsat - vnadir(vsat)
	call vnadir(vsat,vq41)
	vq42(1) = vsat(1)-vq41(1)
	vq42(2) = vsat(2)-vq41(2)
	vq42(3) = vsat(3)-vq41(3)
	vz(1) = vq42(1)
	vz(2) = vq42(2)
	vz(3) = vq42(3)
*
*      tilt = angbtw(vsat, vsat-vz)
	vq41(1) = vsat(1)-vz(1)
	vq41(2) = vsat(2)-vz(2)
	vq41(3) = vsat(3)-vz(3)
	rq41=angbtw(vsat,vq41)
	tilt=rq41
*
	if(tilt.gt.0.)then
*     uopr=vconic(uopr,vz<x>vsat,tilt)
	vq41(1) = vz(2)*vsat(3) - vsat(2)*vz(3)
	vq41(2) = vz(3)*vsat(1) - vsat(3)*vz(1)
	vq41(3) = vz(1)*vsat(2) - vsat(1)*vz(2)
	call vconic(uopr,vq41,tilt,vq42)
	uopr(1) = vq42(1)
	uopr(2) = vq42(2)
	uopr(3) = vq42(3)
*
	end if
*
*
*      view = vcgr - vsat
	vq41(1) = vcgr(1)-vsat(1)
	vq41(2) = vcgr(2)-vsat(2)
	vq41(3) = vcgr(3)-vsat(3)
	view(1) = vq41(1)
	view(2) = vq41(2)
	view(3) = vq41(3)
*
*      sp = view * (uopr <x> vz)
	vq41(1) = uopr(2)*vz(3) - vz(2)*uopr(3)
	vq41(2) = uopr(3)*vz(1) - vz(3)*uopr(1)
	vq41(3) = uopr(1)*vz(2) - vz(1)*uopr(2)
	rq41=view(1)*vq41(1)+view(2)*vq41(2)+view(3)*vq41(3)
	sp=rq41
	 roc = (sp - spl)/(dt-dtl)
	 dtl = dt
	 dt = dtl - sp/roc
	 spl = sp
*      nits = & + 1
	nits=nits+1
*
*      if(nits <= 16) cycle
*
	if(nits.le.16)then
	go to     2
	end if
*
	 dtdiff = dt - dtl
	 go to 900
    2 continue
*
    4 continue
*
*     vslant = vcgr - vsat
	vq41(1) = vcgr(1)-vsat(1)
	vq41(2) = vcgr(2)-vsat(2)
	vq41(3) = vcgr(3)-vsat(3)
	vslant(1) = vq41(1)
	vslant(2) = vq41(2)
	vslant(3) = vq41(3)
*
	dtpass = dt
*
*     gamma = sign( angbtw(vsat,vcgr), -uopr*vslant)
	rq41=uopr(1)*vslant(1)+uopr(2)*vslant(2)+uopr(3)*vslant(3)
	rq42=-rq41
	rq43=sign(angbtw(vsat,vcgr),rq42)
	gamma=rq43
*
	if(abs(gamma).gt.gammax)then
	dtpass = -1.d+0
	go to 102
	end if
*
*     scan = sign(angbtw(vslant,-vz), gamma)
	vq41(1)= -vz(1)
	vq41(2)= -vz(2)
	vq41(3)= -vz(3)
	rq41=sign(angbtw(vslant,vq41),gamma)
	scan=rq41
*
*     slant = rvmag4(vslant)
	rq41=rvmag4(vslant)
	slant=rq41
*
*     vll = vcoord(dt, vsat, 'CLL')
	call vcoord(dt,vsat,'CLL',vq41)
	vll(1) = vq41(1)
	vll(2) = vq41(2)
	vll(3) = vq41(3)
*
	sclat = vll(1)
	sclong = vll(2)
	elev = 90. - abs(gamma) - abs(scan)
	cmode = 'Asc '
*
*     if(veloc*uk < 0.) cmode = 'Desc'
	rq41=veloc(1)*uk(1)+veloc(2)*uk(2)+veloc(3)*uk(3)
*
	if(rq41.lt.0.)then
	cmode='Desc'
	end if
*
	cdt = cd
102   return
*
900   tdiff = 86400.d+0 * dabs(dtdiff)
	write(*, '('' INVNOR; no convergence; err = '', f12.4,
     1  '' sec'')') TDIFF
	call tcivil(dt0, jy,jm,jd,kh,km,sec,nth)
	write(*, '('' csat, dt0: '', a8,8x,i5,4i3,f6.1)')
     x  csat,jy,jm,jd,kh,km,sec
	call tcivil(dtepoc, jy,jm,jd,kh,km,sec,nth)
	write(*, '('' Epoch of Parameters: '', i11,4i3,f6.1/)')
     x  jy,jm,jd,kh,km,sec
	dtpass = -1.d+0
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function ANGBTW(V1,V2)
*
*     To compute the Real*4 angle in degrees between two  Vector*4:
*
	Implicit None
	Double Precision dot,rq81,da1,da2,darg,d0,dtest,dsqrt,dabs
	Double Precision d1,dsign,dc,dacos
	Real angbtw
*
	parameter(dc = 180.d+0/3.141592653589793d+0, d0=0.d+0,
     1  d1=1.d+0)
	double precision vv1(3),vv2(3)
	real  v1(3),v2(3)
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	vv1(1) = v1(1)
	vv1(2) = v1(2)
	vv1(3) = v1(3)
*
	vv2(1) = v2(1)
	vv2(2) = v2(2)
	vv2(3) = v2(3)
*
	dot=vv1(1)*vv2(1)+vv1(2)*vv2(2)+vv1(3)*vv2(3)
	da1=vv1(1)*vv1(1)+vv1(2)*vv1(2)+vv1(3)*vv1(3)
	da2=vv2(1)*vv2(1)+vv2(2)*vv2(2)+vv2(3)*vv2(3)
	darg = da1*da2
*
	if(darg.le.d0)then
	write(*, '('' ANGBTW: an argument is zero '')')
	angbtw = 360.
*
	else
	dtest = dot/dsqrt(darg)
	if(dabs(dtest).gt.d1)dtest=dsign(d1,dtest)
	angbtw = dc*dacos(dtest)
	end if
*
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function NSCAN(kol,card)
*
*     28 Jun 06; corrected error involving minus sign; see below
*     15 Jun 06; modified to accept minus sign
*     Derived from NSCAN.hlf, used in McEDIT
*
*     This function will scan across blanks on an ASCII line
c     to find the next numeric field, from which it will compute
c     a binary value.  kol does not need to be set to the
*     beginning of the next ASCII integer value.
*
	character*1 card,ch
	dimension card(32)
	data n0,n9/48,57/
c     xxxxxxxxxxxxxxxxxxxxxxxxxxx
	mag = 0
	kmag = 0
	jsign = 1
*
	do 10 n = 1,16
	ch = card(kol)
*
*    Find the start of the next integer field:
	if(ch .eq. ' ') then
*     28 Jun 06; changed here:
	if(kmag.eq.0 .and. jsign.gt.0) go to 10
	go to 22
	end if
*
	if(ch .eq. '-') then
	jsign = -1
	go to 10
	end if
*
	ich = ichar(ch)
	if(ich.lt.n0 .or. ich.gt.n9) go to 22
	kmag = 1
	mag = 10*mag + ich - n0
10    kol = kol + 1
*
22    nscan = jsign*mag
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function JULDAY(jyear,month,jday)
*
c     1 Jun 07; logic somewhat modernized
*     27 Jun 98; 2-digit years less than 70 are construed as
c     21st century.   jyear may optionally be 4 digits to avoid
c     ambiguity.
*
c     To get the Julian day number of the day beginning
c     at 12 Z of the given calendar date.  If jyear is given with
c     only two digits, the 20th century is assumed if it is greater
c     than 69, and the 21st century otherwise.  The routine is
c     valid for any date of the Gregorian calendar later than
c     1 Jan 100 AD, and earlier than 1 Jan 4000 AD, since
c     years divisible by 4000 are common years in the Gregorian
c     calendar. If month = 0, then jday is interpreted as the numeric
c     day of the year, 1-366.
cxr   julday
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if(jyear.lt.0 .or. month.lt.0 .or. month.gt.12) go to 160
	if(jday.le.0 .or. jday.gt.366) go to 160
	if(month.ne.0 .and. jday.gt.31) go to 160
*
	jyr = jyear
	if(jyr .lt. 70) jyr = jyr + 2000
*     added 27 Jun 98
	if(jyr .lt. 100) jyr = jyr + 1900
	kxtra = 1
	kwads = jyr/400
	jyrup = 400*kwads
	nquads = (jyr - jyrup)/4
	julian = 1721060 + 146097*kwads + 1461*nquads - nquads/25
c      1721060 is the J.D.N. of 1 Jan 0000 AD.
	if(mod(nquads,25).ne.0 .or. nquads.eq.0) go to 110
c     in a 'short' quadrennium
*
	julian = julian + 1
	kxtra = 0
c     'JULIAN' now points to the first day of the quadrennium con-
c     taining the given date.
*
110   jyrup = jyrup + 4*nquads
	nyears = jyr - jyrup
*
*     changed here 1 Jun 07
	if(nyears.ge.1 .or. (mod(jyr,100).eq.0 .and. mod(jyr,400).ne.0))
     x  then
c      Common year
	 kode = 1115212
	 julian = julian + 365*nyears + kxtra
*
	   else
c      Leap year
	 kode = 1115208
	end if
*
	if(month .le. 1) go to 140
*
	do 130 j = 2,month
	julian = julian + 31 - mod(kode,4)
130   kode = kode/4
*
140   julday = julian + jday - 1
	return
*
160   write(6,161) jyear,month,jday
161   format(' Error in arguments to JULDAY ',3i12)
	julday = 0
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function SINE(xx)
*     Argument in degrees
*
	Implicit Double Precision(d)
	Parameter( dc = 3.141592653589793d+0/180.d+0)
c     xxxxxxxxxxxxxxxxxxxxxxxxxx
	dxx = xx
	ds = dsin(dc*dxx)
	sine = ds
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function COSINE(xx)
*
	Implicit Double Precision(d)
	Parameter( dc = 3.141592653589793d+0/180.d+0)
c     xxxxxxxxxxxxxxxxxxxxxxxxx
	dxx = xx
	ds = dcos(dc*dxx)
	cosine = ds
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function GEODET(gcl)
*
*     To convert geocentric latitude to geodetic:
*
*  Given an elliptical earth with semi-major axis = 1, and
*  semi-minor axis = b, then from the equation of an ellipse,
*  we have
*
*  b2 x2 + y2 = b2                                  (1)
*  y = b sqrt(1 - x2)                               (2)
*  sqrt(1 - x2) = y/b                               (3)
*  dy/dx = -bx/sqrt(1 - x2)        from (2)         (4)
*
*  But the slope dy/dx of a point on the ellipse is the negative
*  reciprocal of the tangent of the geodetic latitude, i.e.
*
*  tan D = sqrt(1 - x2)/bx          where D = geodetic latitude
*  bx tan D = sqrt(1 - x2) = y/b    using (3)
*  b2 tan D = y/x
*
*  On the other hand, y/x is the tangent of the geocentric
c  latitude.  Hence
*
c   b2 tan D = tan C,   or    tan D = (tan C)/b2    (5)
*
c   which can be used to interconvert latitudes in either direction.
*
	Parameter (smajx=6378388., sminx=6356911.946,
     1  b=sminx/smajx, b2 = b*b)
c     These values of semi-major and semi-minor axis were taken
c     from  the Encyclopedia Britannica (see Geodesy).
*
*     These values are also shown as the
*     "International Ellipsoid" as defined in
*     "American Practical Navigator" (Bowditch), Pub. No. 9,
*     Volume I, p. 1119, Defense Mapping Agency.
*
	save
cxr   sine,cosine,arctan
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	gc = abs(gcl)
*
	if(gc.ge.90.)then
	 gd = 90.
	  else
	 tanc = sine(gc)/cosine(gc)
	 gd = arctan(b2, tanc)
	end if
*
	geodet = gd*sign(1., gcl)
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function GEOCEN(gdl)
*
*     See GEODET for explanation of this algorithm.
*
	Implicit None
	Real gdl,gd,abs,gc,tand,sine,cosine,arctan,b2,geocen,sign,
     1  smajx,sminx,b
*
	parameter (smajx=6378388., sminx=6356911.946,
     1  b=sminx/smajx, b2 = b*b)
cxr   sine,arctan,cosine
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	gd = abs(gdl)
*
	if(gd.ge.90.)then
	gc = 90.
	 else
	tand = sine(gd)/cosine(gd)
	gc = arctan(1., b2*tand)
	end if
*
	geocen = gc*sign(1., gdl)
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function ARCTAN(xx,yy)
*
c     This routine was originally written in Meteorological Fortran
c     (MeteFor).  Some of the original code may appear in statements
c     which have been commented out.  The original MeteFor source
c     code may be more self-documenting and readable.
*
c     19 May 06; on laptop
*     SFS; output is in degrees
*     29 Aug 01; see Approximations for Digital Computers, Hastings.
*
	Implicit None
	Double Precision dpi,d1,dpio4,d0,dr2d,dc1,dc3,dc5,dc7,dc9
	Double Precision dc11,dc13,dc15,dx,dy,dabs,dd,dq,da,da2,dr
	Real xx,yy,arctan
*
	Parameter(dpi=3.14159265358979d+0, d1=1.d+0,
     1  dpio4=dpi/4.d+0, d0=0.d+0, dr2d=180.d+0/dpi,
     1  dc1=.9999993329d+0, dc3=-.3332985605d+0,
     2  dc5=.1994653599d+0, dc7=-.1390853351d+0, dc9=.964200441d-1,
     3  dc11=-.559098861d-1, dc13=.218612288d-1, dc15=-.4054058d-2)
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	dx = xx
	dy = yy
*
	if(dabs(dx).le.d0) then
	if(dabs(dy).le.d0) then
	write(*, '('' ARCTAN: both arguments are zero.'')')
	dd = 720.d+0
	end if
*
	dd = 90.d+0
	if(dy.lt.d0)dd = -90.d+0
*
	else
	dq  = dabs(dy/dx)
	da = (dq - d1)/(dq + d1)
	da2 = da*da
	dr = dpio4 + da*(dc1 + da2*(dc3 + da2*(dc5 + da2*(dc7 + da2*(dc9
     1  + da2*(dc11 + da2*(dc13 + da2*dc15)))))))
	dd = dr2d*dr
*
	if(dx.gt.d0)then
	if(dy.lt.d0)dd = -dd
	 else
	dd = 180.d+0 - dd
	if(dy.lt.d0) dd = -dd
	end if
	end if
*
	arctan = dd
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function VERNEQ(dt)
*
cxr   dfrac
*     To return the longitude (positive east) of the Vernal Equinox
c     at the given time dt, expressed either as JDN or MJDN:
*     See page 37  Montenbruck & Pfleger
*
*     Compared with old VERNEQ.for  9 Dec 91;  constant error of
c     .003 deg;  looks OK.
*
	Implicit None
	Double Precision dt,mjd,mjd0,dint,ut,t,gst,gstc,dfrac,ve
	Real verneq
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	mjd = dt - 2400000.5d0
	if(dt .lt. 1721060.d+0) mjd = mjd + 2440588
	mjd0 = dint(mjd)
	ut = 24.d0 * (mjd - mjd0)
	t = (mjd0 - 51544.5d0)/36525.d0
	gst = 6.697374558d0 + 1.0027379093d0 * ut
     1 + (8640184.812866d0+(.093104d0 - 6.2d-6*t)*t)*t/3600.
	gstc = dfrac(gst/24.d0)
	ve = 360. * (1.d0 - gstc)
	if(ve .gt. 180.d0) ve = ve - 360.
	verneq = ve
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function ERAD84(gdlat)
*
*     28 Aug 04; added comment
*     9 Jan 03; uses the spheroid WGS 84
*     5 Jan 03; Liam Gumley reports that MODIS uses the WGS
c     spheroid of 1984, in which the Equatorial radius is
c     6378.137 km, and eccentricity .081819190843, which implies
c     a polar radius of 6356.75231424498 km, i.e.
*     b2 = a2(1 - e2)
*
c     To compute the local radius of the earth in kilometers as a
c     function of geodetic latitude:
*     See EARCOMP.hlf
*
	Implicit Double Precision (d)
*
	Parameter (dmajor=6378.137d+0, dminor=6356.752314245d+0,
     1  b = dminor/dmajor, dmaj2=dmajor*dmajor,
     2  dratio=dminor/dmajor, degrad=3.1415926535898d+0/180.d+0)
*     This is the WGS84 Spheroid.
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if(abs(gdlat) .gt. 90.) then
	write(*, '('' ERAD84: Given Latitude Is Out of Range: '',
     1  1pe16.4)') gdlat
	erad84 = -1.
	end if
*
	gclat = geocen(gdlat)
*     converts to geocentric
	dxe = dmajor * dcos(degrad * gclat)
	dxe2 = dxe*dxe
	dye = dratio*dsqrt(dmaj2 - dxe2)
	erad84 = dsqrt(dxe2 + dye*dye)
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Character*132 Function CAPLET(cvar,length)
c     To convert any lower case ASCII letters to capitals
	character*132 cvar,cout
	character*1 cc,char
	save
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	cout = ' '
	cout(1:length) = cvar(1:length)
*
	do  j = 1,length
	 cc = cout(j:j)
	 ich = ichar(cc)
	 if(ich.ge.97 .and. ich.le.122) cout(j:j) = char(ich-32)
	end do
*
	caplet = cout
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function INPARM(jymd,csat)
*
	implicit double precision(d)
	character*8 csat,cloc
	character*80 cline
	character*132 caplet
	real buf(13),hold(13),de2ra
	logical lexist
	integer jqmax/z'7fffffff'/
	common/e1/dtelem,ds50,xmo,xnodeo,omegao,eo,xincl,xno,xndt2o,xndd60
     x,bstar,x,y,z,xdot,ydot,zdot
	common/c1/ck2,ck4,e6a,qoms2t,s,tothrd,xj3,xke,eqrad,xmnpda,ae
	common/c2/de2ra,pi,pio2,twopi,x3pio2
	equivalence(buf(1),dtelem)
	save
*
	data xj2,xj4,qo,so/1.082616e-3, -1.65597e-6,120., 78./
	data init/0/
cxr   kpattn,caplet,nthday,dtmjdn,isopen
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if(init.eq.0)then
	inquire(file='TLINSET', exist=lexist)
*
	if(.not.lexist)then
	write(*, '('' The Two-Line Element File TLINSET is required,''/
     1  '' but does not exist.''/)')
	dtelem = -1.d+0
	go to 202
	end if
*
	 e6a = 1.e-6
	 pi = 3.14159265
	 de2ra = pi/180.
	 pio2 = pi/2.
	 tothrd = .66666667
	 twopi = 2.*pi
	 x3pio2 = 3.*pio2
	 xj3 = -.253881e-5
	 xke = .743669161e-1
	 eqrad = 6378.135
	 xmnpda = 1440.
	 ae = 1.
	 ck2 = .5*xj2
	 ck4 = -.375*xj4
	 qoms2t = ((qo - so)/eqrad)**4
	 s = (1. + so/eqrad)
	 init = 1
	end if
*
	inp = 1
	cloc = ' '
*
	if(abs(pi-3.14159).gt..001.or.abs(xmnpda-1440.).gt..1)then
	 write(*, '('' INPARM: common blocks c1 and c2 appear not''/
     1  '' to be properly initialized.''/)')
	 write(*,'('' INPARM should be called by VNORAD/VNORAX.'')')
	 go to 202
	end if
*
	do   2 i=1,8
	if(csat(i:i).eq.' ')go to     4
*
	 cloc(i:i) = csat(i:i)
    2 continue
*
    4 continue
	nth = nthday(jymd/10000, mod(jymd/100,100), mod(jymd,100))
	jyyddd = 1000*(jymd/10000) + nth
*
	if(jyyddd.lt.100000)then
	if(jyyddd.gt.70000)jyyddd=jyyddd+1900000
	if(jyyddd.le.70000)jyyddd=jyyddd+2000000
	end if
*
	cloc=caplet(cloc,8)
	if(isopen(58).eq.0)open(58,file='TLINSET',
     x form='formatted',STATUS='old',ACCESS='sequential')
	relmin = 1.e+30
	kfound = 0
	dtymd = dtmjdn(jyyddd/1000, 0, mod(jyyddd,1000), 0,0,0)
*
	do   6 kq0001 = 1,jqmax
	 read(58, '(a80)',end=990,err=992) cline
	 kpat = kpattn(cline,80,cloc,8,0)
	if(kpat.eq.0)go to  6
*
	 kfound = 1
	 read(58,'(18x,d14.8,1x,f10.8,2(1x,f6.5,i2))',end=992,err=992)
     1  dtelem,xndt2o,xndd60,iexp,bstar,ibexp
	 read(58, '(7x, 2(1x,f8.4), 1x, f7.7, 2(1x,f8.4), 1x, f11.8)',
     x  err=992)
     1  xincl,xnodeo,eo,omegao,xmo,xno
	 jepoch = dtelem
	if(jepoch.ge.70000)jepoch=jepoch+1900000
	if(jepoch.lt.70000)jepoch=jepoch+2000000
	 dtep = dtmjdn(jepoch/1000, 0, mod(jepoch,1000),0,0,0)
	 rel = dtymd - dtep
	if(abs(rel).lt.1.)goto20
	if(rel.lt.0..and.abs(rel).gt.relmin)goto990
*
	if(abs(rel).lt.relmin)then
	 relmin = abs(rel)
*
	do   8 i=1,13
8     hold(i) = buf(i)
*
	end if
*
    6 continue
*
	write(*, '('' INPARM: no elements found for '', i10)') jymd
	dtelem = -1.d+0
	go to 202
*
20    xndd60 = xndd60*(10.**iexp)
	xnodeo = xnodeo*de2ra
	omegao = omegao*de2ra
	xmo = xmo*de2ra
	xincl = xincl*de2ra
	temp = twopi/(xmnpda*xmnpda)
	xno = xno*temp*xmnpda
	xndt2o = xndt2o*temp
	xndd60 = xndd60*temp/xmnpda
	a1 = (xke/xno)**tothrd
	temp = 1.5*ck2*(3.*cos(xincl)**2-1.)/(1.-eo*eo)**1.5
	del1 = temp/(a1*a1)
	ao  = a1*(1. - del1*(.5*tothrd + del1*(1. + 134.*del1/81.)))
	delo = temp/(ao*ao)
	xnodp = xno/(1.+delo)
	ideep = 0
	if(twopi/(xnodp*xmnpda).ge..15625)ideep=1
*
	if(ideep.gt.0)then
	write(*, '('' This orbit requires deep space logic'')')
	dtelem = -1.d+0
	go to 202
	end if
*
	bstar = bstar*(10.**ibexp)/ae
	close(58)
	inp = 0
202   inparm = inp
	return
*
990   if(kfound.eq.0)then
	write(*, '('' INPARM: no parameters found for this satellite: '',
     1  a8)') CSAT
	dtelem = -1.d+0
	go to 202
	end if
*
	do  10 i=1,13
10    buf(i) = hold(i)
*
	go to 20
992   write(*, '('' INPARM input error '')')
	kfound = 0
	go to 990
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Subroutine SGP8(iflag,tsince)
*
c     18 Apr 07; NSC
*     8 Feb 99; ARCTAN ===> ACTAN; the result must be augmented
c     by 360 deg, or 2 pi, if it is negative, to conform to the
c     result from ACTAN.
cxr   arctan,fmod2p
*
*     11 Dec 98; received from UW astronomy department, via FTP
c       from sal.wisc.edu; sign in as anonymous; go to the
c       directory pub/jwp; then FTP the file spacetrk.zip;  use
c       the command 'unzip' on Unix
*
*     This routine was written by NORAD, not by SSEC.
*
*     The REAL input argument 'tsince' is time since epoch in
c     minutes.  It may be positive or negative.
*     The epoch itself is the Double Precision variable 'epoch'
c     in the common block /e1/, in the form YYDDD.fffff...,
c     where ffff is a fraction of a day.   The integer 'iflag'
*     should be 1 on the first call to perform
c     needed initialization, and zero thereafter, until and unless
c     another satellite is to be predicted, in which case it must
c     again be 1.
*
*   The outputs from this subroutine are returned thru the
c   common block /e1/.
c   The principal outputs are satellite position (x,y,z),
c   and satellite velocity (xdot,ydot,zdot).
*   The output vectors are celestial, not terrestrial.
*   The position vector uses earth-radius as its unit of length
*   The velocity vector (returned thru common) is in
c   radians/minute.
*   This routine is typically called thru an intermediate subroutine,
c   e.g. like VNORAD.hlf (q.v.).  VNORAD further insures that
c   the needed orbital parameters from the ASCII file TLINSET
c   have been installed in common/e1/.
c     cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	double precision epoch, ds50
	common/e1/epoch,ds50,xmo,xnodeo,omegao,eo,xincl,xno,xndt2o,
     1   xndd6o,bstar,   x,y,z,xdot,ydot,zdot
	common/c1/ck2,ck4,e6a,qoms2t,s,tothrd,
     1   xj3,xke,xkmper,xmnpda,ae
	save
*
	data rho/.15696615/
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if (iflag .eq. 0) go to 100
*
*    Recover original mean motion (xnodp) and semimajor axis (aodp)
*    from input elements --------- calculate ballistic coefficient
*    (b term) from input b* drag term
*
	a1 = (xke/xno)**tothrd
*     xno is mean motion in radians/minute
	cosi = cos(xincl)
	theta2 = cosi*cosi
	tthmun = 3.*theta2 - 1.
	eosq = eo*eo
	betao2 = 1.-eosq
	betao = sqrt(betao2)
	del1=1.5*ck2*tthmun/(a1*a1*betao*betao2)
*     ao=a1*(1.-del1*(.5*tothrd+del1*(1.+134./81.*del1)))
	ao = a1*(1. - del1*(.5*tothrd + del1*(1. + 134.*del1/81.)))
	delo = 1.5*ck2*tthmun/(ao*ao*betao*betao2)
	aodp = ao/(1.-delo)
	xnodp = xno/(1.+delo)
	b = 2.*bstar/rho
*
*      Initialization
*
	isimp = 0
	po = aodp*betao2
	pom2 = 1./(po*po)
	sini = sin(xincl)
	sing = sin(omegao)
	cosg = cos(omegao)
	temp = .5*xincl
	sinio2 = sin(temp)
	cosio2 = cos(temp)
	theta4 = theta2**2
	unm5th = 1.-5.*theta2
	unmth2 = 1.-theta2
*     a3cof = -xj3/ck2*ae**3
	a3cof = -xj3 * ae**3 / ck2
	pardt1 = 3.*ck2*pom2*xnodp
	pardt2 = pardt1*ck2*pom2
	pardt4 = 1.25*ck4*pom2*pom2*xnodp
	xmdt1 = .5*pardt1*betao*tthmun
	xgdt1 = -.5*pardt1*unm5th
	xhdt1 = -pardt1*cosi
	xlldot = xnodp + xmdt1 +
     2   .0625*pardt2*betao*(13.-78.*theta2+137.*theta4)
	omgdt=xgdt1+
     1  .0625*pardt2*(7.-114.*theta2 + 395.*theta4) + pardt4*(3.-36.*
     2   theta2+49.*theta4)
	xnodot=xhdt1+
     1   (.5*pardt2*(4.-19.*theta2)+2.*pardt4*(3.-7.*theta2))*cosi
	tsi = 1./(po-s)
	eta = eo*s*tsi
	eta2 = eta**2
	psim2 = abs(1./(1.-eta2))
	alpha2 = 1.+eosq
	eeta = eo*eta
	cos2g = 2.*cosg**2-1.
	d5 = tsi*psim2
	d1 = d5/po
	d2 = 12. + eta2*(36. + 4.5*eta2)
	d3 = eta2*(15. + 2.5*eta2)
	d4 = eta*(5. + 3.75*eta2)
	b1 = ck2*tthmun
	b2 = -ck2*unmth2
	b3 = a3cof*sini
	c0 = .5*b*rho*qoms2t*xnodp*aodp*tsi**4*psim2**3.5/sqrt(alpha2)
	c1 = 1.5*xnodp*alpha2**2*c0
	c4 = d1*d3*b2
	c5 = d5*d4*b3
	xndt=c1*(
     1  (2. + eta2*(3. + 34.*eosq) + 5.*eeta*(4. + eta2) + 8.5*eosq) +
     1   d1*d2*b1 +  c4*cos2g+c5*sing)
	xndtn=xndt/xnodp
*
*      If drag is very small, the isimp flag is set and the
*      equations are truncated to linear variation in mean
*      motion and quadratic variation in mean anomaly
	if(abs(xndtn*xmnpda) .lt. 2.16e-3) go to 50
*
	d6 = eta*(30. + 22.5*eta2)
	d7 = eta*(5. + 12.5*eta2)
	d8 = 1. + eta2*(6.75 + eta2)
	c8 = d1*d7*b2
	c9 = d5*d8*b3
	edot = -c0*(
     1   eta*(4.+eta2+eosq*(15.5+7.*eta2))+eo*(5.+15.*eta2)+
     1   d1*d6*b1 +
     1   c8*cos2g+c9*sing)
	d20 = .5*tothrd*xndtn
	aldtal = eo*edot/alpha2
	tsdtts = 2.*aodp*tsi*(d20*betao2 + eo*edot)
	etdt = (edot + eo*tsdtts)*tsi*s
	psdtps = -eta*etdt*psim2
	sin2g = 2.*sing*cosg
	c0dtc0 = d20+4.*tsdtts-aldtal-7.*psdtps
	c1dtc1 = xndtn + 4.*aldtal + c0dtc0
	d9 = eta*(6. + 68.*eosq) + eo*(20. + 15.*eta2)
	d10 = 5.*eta*(4. + eta2) + eo*(17. + 68.*eta2)
	d11 = eta*(72. + 18.*eta2)
	d12 = eta*(30. + 10.*eta2)
	d13 = 5. + 11.25*eta2
	d14 = tsdtts-2.*psdtps
	d15 = 2.*(d20 + eo*edot/betao2)
	d1dt = d1*(d14+d15)
	d2dt = etdt*d11
	d3dt = etdt*d12
	d4dt = etdt*d13
	d5dt = d5*d14
	c4dt = b2*(d1dt*d3+d1*d3dt)
	c5dt = b3*(d5dt*d4+d5*d4dt)
	d16=
     1     d9*etdt+d10*edot +
     1     b1*(d1dt*d2+d1*d2dt) +
     1     c4dt*cos2g+c5dt*sing+xgdt1*(c5*cosg-2.*c4*sin2g)
	xnddt = c1dtc1*xndt + c1*d16
	eddot = c0dtc0*edot-c0*(
     1   (4.+3.*eta2+30.*eeta+eosq*(15.5+21.*eta2))*etdt+(5.+15.*eta2
     x    +eeta*(31.+14.*eta2))*edot +
     1   b1*(d1dt*d6 + d1*etdt*(30. + 67.5*eta2))   +
     1   b2*(d1dt*d7 + d1*etdt*(5. + 37.5*eta2))*cos2g +
     1   b3*(d5dt*d8 + d5*etdt*eta*(13.5 + 4.*eta2))*sing + xgdt1*(c9*
     x      cosg-2.*c8*sin2g))
	d25 = edot**2
	d17 = xnddt/xnodp - xndtn**2
	tsddts=2.*tsdtts*(tsdtts-d20)+aodp*tsi*(tothrd*betao2*d17-4.*d20*
     x   eo*edot+2.*(d25+eo*eddot))
	etddt =(eddot+2.*edot*tsdtts)*tsi*s+tsddts*eta
	d18 = tsddts-tsdtts**2
	d19= - psdtps**2/eta2 - eta*etddt*psim2 - psdtps**2
	d23 = etdt*etdt
	d1ddt=d1dt*(d14+d15)+d1*(d18-2.*d19+tothrd*d17+2.*(alpha2*d25
     x   /betao2+eo*eddot)/betao2)
	xntrdt=xndt*(2.*tothrd*d17+3.*
     1  (d25+eo*eddot)/alpha2 - 6.*aldtal**2 +
     1  4.*d18-7.*d19 ) +
     1  c1dtc1*xnddt+c1*(c1dtc1*d16+
     1  d9*etddt+d10*eddot+d23*(6.+30.*eeta+68.*eosq)+
     1  etdt*edot*(40.+30.*
     x  eta2+272.*eeta)+d25*(17.+68.*eta2) +
     1    b1*(d1ddt*d2+2.*d1dt*d2dt+d1*(etddt*d11+d23*(72.+54.*eta2))) +
     1    b2*(d1ddt*d3+2.*d1dt*d3dt+d1*(etddt*d12+d23*(30.+30.*eta2))) *
     1    cos2g+
     1      b3*((d5dt*d14+d5*(d18-2.*d19)) *
     1 d4+2.*d4dt*d5dt+d5*(etddt*d13+22.5*eta*d23)) *sing+xgdt1*
     1       ((7.*d20+4.*eo*edot/betao2)*
     x    (c5*cosg-2.*c4*sin2g)
     x      + ((2.*c5dt*cosg-4.*c4dt*sin2g)-xgdt1*(c5*sing+4.*
     x      c4*cos2g))))
	tmnddt = xnddt*1.e9
	temp = tmnddt**2-xndt*1.e18*xntrdt
	pp = (temp+tmnddt**2)/temp
	gamma  =  -xntrdt/(xnddt*(pp-2.))
	xnd = xndt/(pp*gamma)
	qq = 1.-eddot/(edot*gamma)
	ed = edot/(qq*gamma)
	ovgpp = 1./(gamma*(pp+1.))
	go to 70
*
   50 isimp=1
	edot=-tothrd*xndtn*(1.-eo)
*
   70 iflag=0
*
*      Update for secular gravity and atmospheric drag
*
  100 xmam = fmod2p(xmo + xlldot*tsince)
	omgasm = omegao + omgdt*tsince
	xnodes = xnodeo + xnodot*tsince
	if(isimp .eq. 1) go to 105
*
	temp = 1.-gamma*tsince
	temp1 = temp**pp
	xn = xnodp+xnd*(1.-temp1)
	em = eo+ed*(1.-temp**qq)
	z1 = xnd*(tsince+ovgpp*(temp*temp1-1.))
	go to 108
*
  105 xn = xnodp+xndt*tsince
	em = eo+edot*tsince
	z1 = .5*xndt*tsince*tsince
*
  108 z7 = 3.5*tothrd*z1/xnodp
	xmam = fmod2p(xmam+z1+z7*xmdt1)
	omgasm = omgasm+z7*xgdt1
	xnodes = xnodes+z7*xhdt1
*
*      Solve Keplers equation
	zc2=xmam+em*sin(xmam)*(1.+em*cos(xmam))
*
	do 130 i=1,10
	sine=sin(zc2)
	cose=cos(zc2)
	zc5=1./(1.-em*cose)
	cape=(xmam+em*sine-zc2)*
     1   zc5+zc2
	if(abs(cape-zc2) .le. e6a) go to 140
*
  130 zc2=cape
*
*      short period preliminary quantities
*
  140 am = (xke/xn)**tothrd
	beta2m = 1.-em*em
	sinos = sin(omgasm)
	cosos = cos(omgasm)
	axnm = em*cosos
	aynm = em*sinos
	pm = am*beta2m
	g1 = 1./pm
	g2 = .5*ck2*g1
	g3 = g2*g1
	beta = sqrt(beta2m)
	g4 = .25*a3cof*sini
	g5 = .25*a3cof*g1
	snf = beta*sine*zc5
	csf = (cose-em)*zc5
*     fm = actan(snf,csf)
	fm = .0174532925*arctan(csf,snf)
	if(fm .lt. 0.) fm = fm + 6.2831853
*     To simulate ACTAN, which returns 0 to 2*pi
	snfg = snf*cosos + csf*sinos
	csfg = csf*cosos - snf*sinos
	sn2f2g = 2.*snfg*csfg
	cs2f2g = 2.*csfg**2-1.
	ecosf = em*csf
	g10 = fm - xmam+em*snf
	rm = pm/(1.+ecosf)
	aovr = am/rm
	g13 = xn*aovr
	g14 = -g13*aovr
	dr = g2*(unmth2*cs2f2g - 3.*tthmun) - g4*snfg
	diwc = 3.*g3*sini*cs2f2g-g5*aynm
	di = diwc*cosi
*
*      update for short period periodics
*
	sni2du = sinio2*(
     1   g3*(.5*(1.-7.*theta2)*sn2f2g-3.*unm5th*g10)-g5*sini*csfg*(2.+
     2      ecosf)) - .5*g5*theta2*axnm/cosio2
	xlamb=fm+omgasm+xnodes+g3*(.5*(1.+6.*cosi-7.*theta2)*sn2f2g-3.*
     1   (unm5th+2.*cosi)*g10)+g5*sini*(cosi*axnm/(1.+cosi)-(2.
     2    + ecosf)*csfg)
	y4 = sinio2*snfg+csfg*sni2du+.5*snfg*cosio2*di
	y5 = sinio2*csfg-snfg*sni2du+.5*csfg*cosio2*di
	r = rm + dr
	rdot = xn*am*em*snf/beta + g14*(2.*g2*unmth2*sn2f2g+g4*csfg)
	rvdot = xn*am**2*beta/rm +
     1    g14*dr+am*g13*sini*diwc
*
*      orientation vectors
*
	snlamb = sin(xlamb)
	cslamb = cos(xlamb)
	temp = 2.*(y5*snlamb-y4*cslamb)
	ux = y4*temp + cslamb
	vx = y5*temp-snlamb
	temp = 2.*(y5*cslamb + y4*snlamb)
	uy = -y4*temp + snlamb
	vy = -y5*temp + cslamb
	temp = 2.*sqrt(1. - y4*y4 - y5*y5)
	uz = y4*temp
	vz = y5*temp
*
*     Position and Velocity
	x = r*ux
	y = r*uy
	z = r*uz
	xdot = rdot*ux + rvdot*vx
	ydot = rdot*uy + rvdot*vy
	zdot = rdot*uz + rvdot*vz
*     Note: x,y,z are in units of earth radius.
c     xdot,ydot,zdot are in radians per minute   - Nagle
*
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function SMAXIS(xno,ep,fincl)
*
*     2 Apr 07; NSC
*     5 Oct 99; improved accuracy in computing sio
c     Semi-major axis from mean motion, inclination & eccentricity
c     The output is in nautical miles.
*     The inverse function is DMMOTN.hlf
*
c     xno - mean motion in revolutions per day
c     ep - eccentricity
c     fincl - inclination
*     The International Nautical Mile is 1852 meters, exactly.
*     This is 6076.11549 feet, approximately.
*     See Bowditch, Vol II, page 446.
*
	implicit double precision (a-h,o-z)
*
c     Obtained from US Air Force, 1000th Satellite Operations Group,
c     Offutt AFB, Nebraska   68113, Tel 402-294-5455
c     (Attn: Capt Wiseman  -  27 Jul 87)
*
	parameter( const = -1.5d+0 * 1.082616d-3,
     1  rc2 = 20925639.76d+0 ** 2 , third = 1.d+0/3.d+0,
     2  frac = 44.d+0/81.d+0 , ft2nm = 1.64578833d-4, d1=1.d+0)
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	xnonew = 6.2831853071796d+0 * xno / 86400.d+0
	sio = dsin(.0174532925199433d+0 * fincl)
c     1.74532925199433D-02  = pi/180
	ao = ( 1.40764544d+16/xnonew ** 2) ** third
	do =  const * (rc2/ao**2)*(d1 - ep**2) **(-1.5d+0)*
     1  (d1 - 1.5d+0 * sio**2)
	smaxis =
     1  ao * (d1 + do*(third - do*(third - do*frac))) * ft2nm
*     Note: This should probably be multiplied by .999408
*     based on a McIdas test 11 Oct 89.
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
c     Vector Function UNIT4(VEC)
	Subroutine Unit4(vec,uvec)
*
*     2 Apr 03; to compute the norm of an arbitrary vector:
*
	real vec(3),uvec(3)
	equivalence(jr2,r2),(sum2,jsum2)
*     data j$01/z'20000000'/             dec 536870912
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	sum2 = vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)
*
	if(sum2 .gt. 0.)then
	jr2 = jsum2/2 + 536870912
*
	do 2 n = 1,6
	r1 = r2
2     r2 = .5*(r1 + sum2/r1)
*
	do 4 j = 1,3
4     uvec(j) = vec(j)/r2
*
	    else
	do 6 j = 1,3
6     uvec(j) = 0.
*
	write(*, '('' UNIT4: Warning - the given vector argument has '',
     1  ''zero length.'')')
	end if
*
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Subroutine VCONIC(vg,vplane,rot,vqonic)
*
c     This routine was originally written in Meteorological Fortran
c     (MeteFor), an extension of Fortran-77, in order to utilize the
c     vector and matrix notation available in MeteFor.  This
c     routine is also maintained in MeteFor.  Some of the original
c     MeteFor code may appear in statements which have been commented
c     out.  The original MeteFor source (suffix .hlf) is more
c     readable and self-documenting.  See the document MeteFor.doc.
*
*     21 Aug 98; shortened slightly
*     Imagine a cone whose axis is the vector VCONE.  Let the
c     vector VG be an element of the cone directed away from the
c     apex of the cone.  This vector-valued function rotates the
c     vector VG rightward (clockwise) about the axis VCONE, as
c     VCONE is pointed toward the viewer. If the cone is degenerate,
c     i.e. if it is a plane, then this function is the same as
c     VROT4.
*
c       More crudely, imagine looking into the interior of an empty
c     ice cream cone.  The vector VCONE is the axis of the cone,
c     directed toward the viewer.  Let the vector VG
c     be any element of the cone extending from the apex of
c     the cone toward the lip.  If the viewer now rotates the cone
c     'rot' degrees clockwise from his perspective, this function
c     computes the new vector value of VG after the rotation.
*     If rot<0, then the rotation is counter-clockwise.
*     The cone may be merely a plane, in which case VCONE and VG are
c     normal.  A practical application would be determining the direction
c     of view of a conically-scanning satellite.
*
	Implicit None
	real vq41(3),vq42(3),vq43(3),vq44(3),vq45(3)
	real  uv(3),vplane(3),ur(3),vg(3),uh(3),sine,cosine,rot
	real rq41,rq42,rq43,rq44,vqonic(3)
cxr   sine,cosine,unit4
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
*     uv = unit4(vplane)
	call unit4(vplane,uv)
*
*     ur = unit4(vg ** uv)
	vq41(1) = vg(2)*uv(3) - uv(2)*vg(3)
	vq41(2) = vg(3)*uv(1) - uv(3)*vg(1)
	vq41(3) = vg(1)*uv(2) - uv(1)*vg(2)
	call unit4(vq41,ur)
*
*     uh = unit4(uv ** ur)
	vq41(1) = uv(2)*ur(3) - ur(2)*uv(3)
	vq41(2) = uv(3)*ur(1) - ur(3)*uv(1)
	vq41(3) = uv(1)*ur(2) - ur(1)*uv(2)
	call unit4(vq41,uh)
*
*     vconic = (uh*vg) * (ur*sine(rot) + uh*cosine(rot)) +
*    1  (vg*uv) * uv
	rq41=sine(rot)
	rq42=cosine(rot)
	rq43=uh(1)*vg(1)+uh(2)*vg(2)+uh(3)*vg(3)
	vq41(1) = ur(1)*rq41
	vq41(2) = ur(2)*rq41
	vq41(3) = ur(3)*rq41
	vq42(1) = uh(1)*rq42
	vq42(2) = uh(2)*rq42
	vq42(3) = uh(3)*rq42
	vq43(1) = vq41(1)+vq42(1)
	vq43(2) = vq41(2)+vq42(2)
	vq43(3) = vq41(3)+vq42(3)
	rq44 = vg(1)*uv(1)+vg(2)*uv(2)+vg(3)*uv(3)
	vq44(1) = vq43(1)*rq43
	vq44(2) = vq43(2)*rq43
	vq44(3) = vq43(3)*rq43
	vq45(1) = uv(1)*rq44
	vq45(2) = uv(2)*rq44
	vq45(3) = uv(3)*rq44
	vqonic(1) = vq44(1)+vq45(1)
	vqonic(2) = vq44(2)+vq45(2)
	vqonic(3) = vq44(3)+vq45(3)
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function ARCCOS(xx)
*
*     19 Feb 03; tested superfically in NORWATCH
	implicit double precision(d)
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	x = xx
*
	if(abs(x).gt.1.)then
	write(*, '('' ARCCOS: argument > 1: '', 1pe14.6)') c
	ac = 999.
*
	else
	dx = x
	da = 57.2957795131*dacos(dx)
	ac = da
	end if
*
	arccos = ac
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function EARRAD(gdlat)
*
*     9 Jan 03; uses the International Spheroid (see Bowditch p. 1194)
*     10 Aug 01; geodetic ===> geocentric as input
c     To compute the local radius of the earth in kilometers as a
c     function of geodetic latitude:
*
	Implicit Double Precision (d)
	Parameter (dmajor=6378.388d+0, dminor=6356.911946d+0,
     1  b = dminor/dmajor, dmaj2=dmajor*dmajor,
     2  dratio=dminor/dmajor, degrad=3.1415926535898d+0/180.d+0)
*     Semi-major and -minor axes of Earth in km
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if(abs(gdlat) .gt. 90.) then
	write(*, '('' EARRAD: Given Latitude Is Out of Range: '',
     1  1pe16.4)') gdlat
	earrad = -1.
	end if
*
	gclat = geocen(gdlat)
*     converts to geocentric
	dxe = dmajor * dcos(degrad * gclat)
	dxe2 = dxe*dxe
	dye = dratio*dsqrt(dmaj2 - dxe2)
	earrad = dsqrt(dxe2 + dye*dye)
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function RVMAG4(VEC)
*     to compute the length of a Vector*4:
*
	real vec(3)
	equivalence(sum,jsum), (r1,jr1)
c     xxxxxxxxxxxxxxxxxxxxxxxxxx
	r1 = 0.
	sum = vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)
*
	if(sum.gt.0.) then
	jr1 = jsum/2 + 536870912
*
	do 2 n = 1,6
	r0 = r1
2     r1 = .5*(sum/r0 + r0)
*
	end if
*
	rvmag4 = r1
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Double Precision Function DFRAC(ditem)
*
*     The function returns the fraction from the algebraically
*     lesser integer, i.e.
*
c     1.7 ===> .7   and   -1.7 ===> .3
*
	Implicit Double Precision (a-h,o-z)
cxr   dflt
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	df = ditem
*
	if(dabs(df) .gt. 2147483647.9999d0) then
	write(*, '('' DFRAC: argument is too large: '', d20.10)') df
	dfrac = 999999999.
*
	  else
	item = df
	df = df - dflt(item)
	if(df .lt. 0.d+0) df = df + 1.d+0
	dfrac = df
	end if
*
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function KPATTN(cline,lenli,cpat,leng,kode)
*
*     To find the column (1-lenli) where a given pattern begins
*     If kode=0, searches for exact pattern in 'cpat'
c     If kode!=0, disregards lower/upper case, i.e.  xyz = XYZ
*     If the pattern does not exist, zero is returned.
*     leng  is the length of the pattern, with maximum 24.
*
	character*132 cline,caplet
	character*24 csub,cpatl,cpat
cxr   caplet
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	len = leng
	if(len.gt.24)len=24
	cpatl = ' '
	cpatl = cpat(1:len)
	if(kode.ne.0)cpatl=caplet(cpatl,len)
*
	do   2 k=1,lenli
	 k1 = k
	 k2 = k1 + len - 1
*
	if(k2.gt.lenli)go to 4
*
	 csub = cline(k1:k2)
	if(kode.eq.1)csub=caplet(csub,len)
	if(csub.eq.cpatl) go to 20
*
    2 continue
*
    4 continue
	kpattn = 0
	return
*
20    kpattn = k1
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function NTHDAY(ny,nm,nd)
*
c     31 May 04; if,then  ===> go to
*     This function returns the day of the year (1-366), given
c     the civil year,month,day.
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	kode = 1115212
	if(mod(ny,4).eq.0 .and. mod(ny,100).ne.0) kode = 1115208
	if(mod(ny,400) .eq. 0) kode = 1115208
	nth = nd
*
	if(nm .gt. 1) then
	nth = 0
	m = nm - 1
*
	do 110 j = 1,m
	nth = nth + 31 - mod(kode,4)
110   kode = kode/4
*
	nth = nth + nd
	end if
*
	nthday = nth
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function ISOPEN(lun)
*     To ascertain if logical unit number lun is open:
c     (1=yes, 0=no)
*
	logical lopen
c     xxxxxxxxxxxxxxxxxxxxx
	inquire(lun, opened=lopen)
	iset = 0
	if(lopen) iset = 1
	isopen = iset
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Function FMOD2P(x)
*     Used in connection with NORAD SGP and SGP4
*     This routine was not written by SSEC
*
	common/c2/de2ra,pi,pio2,twopi,x3pio2
	save
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	fmod2p = x
	i = fmod2p/twopi
	fmod2p = fmod2p-i*twopi
	if(fmod2p.lt.0) fmod2p = fmod2p+twopi
	return
	end
*      rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Double Precision Function DFLT(narg)
c     DFLOAT does not work on some Unix platforms
	Double Precision dd
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	dd = narg
	dflt = dd
	return
	end
c     rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Subroutine TCIVIL(dtime,  nyr,nmo,nda,nhr,min,sec,  nthday)
*
*     SFS
*     TCIVIL - to convert time in JDN or MJDN to time
c     expressed in civil units.  Nthday varies from 1 to 366.
*     This differs from TINVER in that seconds are returned as
c     REAL*4, and hence may be fractional.
c     A four-digit year nyr is returned.
*
	Implicit None
	Double Precision dtime,dt,dfrac,dflt
	Real secs,fksecs,frac,fnsec,sec
	Integer jdt,jdtx,nyr,nmo,nda,nthday,ksecs,nhr,min,nsec
*
	character*4 cmon
cxr   julinv,dflt
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	if(dtime .le. 0.d+0) go to 800
*
	dt = dtime + .5d+0
	jdt = dt
	jdtx = jdt
	if(dt.lt.1721060.d+0)jdtx = jdtx + 2440588
*     An input smaller than 1721060 (beginning of Christian era)
c     is presumed to be a Modified Julian Day Number.
	call julinv(jdtx,  nyr,nmo,nda,  cmon,nthday)
	dfrac = dt - dflt(jdt)
	secs = 86400.d+0 * dfrac + 1.d-6
	ksecs = secs
	nhr = ksecs/3600
	min = (ksecs - 3600*nhr)/60
	nsec = ksecs - 3600*nhr - 60*min
	fksecs = ksecs
	frac = secs - fksecs
	fnsec = nsec
	sec = fnsec + frac
	if(sec.ge.60.)sec=59.99999
	return
*
800   nyr = 0
	nmo = 0
	nda = 0
	nhr = 0
	min = 0
	sec = 0
	write(*,'('' WARNING: Input to TCIVIL is bad: '',1pd20.10)')dtime
	return
	end
c     rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
	Subroutine JULINV(julin, ky,km,kd, cmon,nth)
*
*     To convert the Julian Day Number (JDN) to civil date:
*     7 Aug 98; tested using JULTEST
*     This routine works for years 4 AD to 3999 AD.
*     A 4-digit year ky is returned.  The returned ASCII month is of
*     the form ' Jan', i.e. character*4 with a lead blank.  nth
*     is the day of the year (1-366).
*
	Implicit None
	Integer jdn,julin,ndce,kwads,nrem,long,ncents,kq0001,nquads
	Integer kq0002,kyear,kq0003,kym1,leap,mod,kode,nth,mm,m
	Integer ndm,kd,km,ky
	character*4 cmon,cmons(12)
	integer jqmax/z'7fffffff'/
*
	data cmons/' Jan',' Feb',' Mar',' Apr',' May',' Jun',' Jul',
     1  ' Aug',' Sep',' Oct',' Nov',' Dec'/
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	jdn = julin
	if(jdn<1721060 .or. jdn>3182410) go to 900
*
*     from beginning of Christian era to about 4000 AD
	ndce = jdn - 1721060
	kwads = ndce/146097
	nrem = ndce - 146097*kwads
	long = 1
	ncents = 0
*
	do   2 kq0001 = 1,jqmax
*      nrem = & - (36524 + long)
	 nrem = nrem-(36524+long)
	 if(nrem.lt.0) go to   4
*
	 long = 0
2     ncents = ncents+1
*
4     nrem = nrem + (36524 + long)
	nquads = 0
*
	do   6 kq0002 = 1,jqmax
	nrem = nrem-(1460+long)
	if(nrem.lt.0) go to   8
*
	 long = 1
6     nquads = nquads+1
*
8     nrem = nrem + (1460 + long)
	kyear = 400*kwads + 100*ncents + 4*nquads
	long = 1
*
	do  10 kq0003 = 1,jqmax
	nrem = nrem-(365+long)
	if(nrem.lt.0) go to  12
*
	 long = 0
10    kyear = kyear+1
*
12    nrem = nrem+(365+long)
	kym1 = kyear - 1
	leap = (mod(kym1,4)+1)/4 - (mod(kym1,100)+1)/100
     1  + (mod(kym1,400)+1)/400
	kode = 1115212 - 4*leap
	nth = nrem + 1
*
	do  14 m = 1,12
	 mm = m
	 ndm = 31 - mod(kode,4)
	 nrem = nrem-ndm
	 if(nrem.lt.0) go to  16
*
14    kode = kode/4
*
16    nrem = nrem+ndm
	kd = nrem + 1
	km = mm
	cmon = cmons(mm)
	ky = kyear
	return
*
900   write(*, '('' JULINV: JDN out of range:'', i14)') jdn
	return
	end