# Routine Descriptions

## ABERANNUAL - COMPUTE THE ANNUAL ABERRATION OFFSET.

[Next Routine] [List of Routines]
```NAME:
aberAnnual - compute the annual aberration offset.
SYNTAX: v=aberAnnual(julianDate)
ARGS:
juliandate[n]: double julian data of interest. ut1 or utc base is equiv.
RETURNS:
v[3,n]: double 3vector for the aberation correction. Add this
to the curret true ra,dec when going ra,dec->az,el.

DESCRIPTION

Return the annual aberration vector offset. This vector should be added
to the true object position to get the apparent position of the object.

Aberration requires the observer to point at a direction different from
the geometrical position of the object. It is caused by the relative motion
of the two objects (think of it as having to bend the umbrella toward your
direction of motion to not get wet..). It is the velocity vector offset of
the observer/object from the geometric direction measured in an inertial
frame.

Complete aberration is called planetary aberration and includes the motion
of the object and the motion of the observer. For distant objects the
objects motion can be ignored. Stellar aberration is the portion due to the
motion of the observer (relative to some inertial frame). It consists of
secular aberration (motion of the solar system), annual aberration
(motion of the earth about the sun), and diurnal aberration (caused by
the rotation of the earth about its axis.

Secular aberration is small and ignored (the solar system is assumed to be
in  the inertial frame of the stars). The earths velocity about the sun of
+/- 30km/sec = +/- .0001v/c gives  +/- 20" variation for annual aberration.
diurnal variation gives a maximum of about +/- .32" at the equator over a day.
This routine computes annual aberration using the approximate routines no
page B17 and C24 of the AA.

EXAMPLE:
```

(See /pkg/rsi/local/libao/phil/gen/pnt/aberannual.pro)

## AGCAZZADIF - COMPUTE AZ,ZA OFFSETS FROM AGC DATA.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
agcazzadif - compute az,za offsets from agc data.

ARGS:
rcvnum  :int      reciever number
raJ     :double   J2000 ra in hours
decJ    :double   J2000 dec in hours

agcDat[4,n]:fltarr   holds tmsecs,azPos,grPos,chPos in degrees
stYear: long    year for starting time to use
stdayNum: long    daynumber of year for starting time to use
stSec: double  start sec
tmStepSec: double  time step we want
npts: long    number of points to compute
KEYWORDS:
raDecPacked: if set the ra is hhmmss.s and dec is ddmmss.s

RETURNS:
azdif[npts]:double offset in azimuth (great circle arc minutes)
zadif[npts]:double offset in za (great circle arc minutes)
stat       : int   1 ok, 0 no data

DESCRIPTION:
When doing on the fly mapping, the telescope is moving rapidly while the
correlator, ri is recording the data. The pointing system can be requested
to log the az,za positions in a file at a 25 hz rate. This routine will
then compute the actual great circle az,za offset of the correlator
sampled data from the center of the map. To do this the user must input:

1. The center of the map in ra,dec J2000.. raJ, decJ.
2. The 25 hz sampled pointing data ... agcDat.
3. The start time of each correlator dump. You enter the
starting year,dayno, astSecond of the scan: stYear,stDay,stSec
as well as the integration time for each record (solar secs): tmStepSec
Finaly you must enter the number of records: npts.

The routine converts raJ,decJ to az,za for the center of each integration
using ao_radecjtoazza. It then interpolates the measured az,za positions
to these times and computes the az,za differences.

Each sample of the 25Hz agc data is stored as 4 longs:
agc.tm  long milliseconds from midnite (AST)
agc.az  long azimuth position in degrees  *10000 (dome side)
agc.gr  long za position of dome in degrees  *10000
agc.ch  long za position of ch in degrees  *10000

To input all the data of a file:
filename='/share/obs1/pnt/pos/dlogPos.1'
openr,lun,filename,/get_lun
fstat=fstat(lun)
numSamples=fstat.size/(4*4)
allocate array to hold all the data
inparr=lonarr(4,numSamples)
free_lun,lun
convert to seconds and  degrees.
datArr=fltarr(4,numSamples)
datArr[0,*]= inpArr[0,*]*.001   ; millisecs to secs
datArr[1,*]= inpArr[1,*]*.0001  ; az data to deg
datArr[2,*]= inpArr[2,*]*.0001  ; gr data to deg
datArr[3,*]= inpArr[3,*]*.0001  ; ch data to deg

NOTE:
The data in agcDat is stored with only a secsFrom midnite timestamp.
The first entry in agcDat should be the same day as stDay. If that
is true,then it is ok to cross midnite.
```

## ALFABMPOS - COMPUTE ALFA BEAM POSITIONS FROM AZ,ZA,JD BEAM 0.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
alfabmpos - compute alfa beam positions from az,za,jd beam 0.
SYNTAX: alfabmpos,az,za,juldat,rahrs,decDeg,rotAngle=rotAngle,\$
hornOffsets=hornOffsets,offsetsonly=offsetsonly,\$
nomodel=nomodel
ARGS:
az[n]:  float  azimuth of central pixel in degrees
za[n]:  float  zenith angle of central pixel in degrees
juldat[n]: double julian date for each az,za. Should include the fraction
of the day.
KEYWORDS:
rotAngl: float  rotation angle (deg)of alfa rotator. default is 0 degrees.
sitting on the rotator floor, positive is clockwise.
offsetsonly:       if set, then only return the offsets, don't bother to
compute the ra,decs
nomodel   :        If set , then then az,za supplied already has the
model removed. It will not remove it a second time.
RETURNS:
raHrs[7,n]: float ra in hours for the 7  beams and the n positions.
decDeg[7,n]: float declination in degrees for the 7  beams and the
n positions.
hornOffsets[2,7]:float The offsets of the horns relative to pixel 0 in
great circle degrees. The first dimension is
az,za. The second dimension is pixel number

DESCRIPTION:
Given the az, za, and juldate for the center pixel of alfa, this routine
will compute the ra,dec for each of the 7 beams of alfa. The order returned
is pixel0,1,2,3,4,5,6,7. If az,za,juldat are arrays of length n then
the return data will be [7,n].
By default the rotation angle for the array is 0 deg. You can uset a
different rotation angle using the rotangle keyword.
The horn offsets relative to pixel 0 are returned in the hornOffsets
keyword. The units are arcseconds. They are the azOff,zaOff values that are

```

(See /pkg/rsi/local/libao/phil/gen/pnt/alfabmpos.pro)

## ANGLESTOVEC3 - CONVERT FROM ANGLES TO 3 VECTORS.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
anglestovec3 - convert from angles to 3 vectors.

SYNTAX: v=anglestovec3(c1Rd,c2Rd)

ARGS:
c1Rd[npts] : float/double first  angular coordinate in radians.
c2Rd[npts] : float/double second angular coordinate in radians.

RETURNS
v[3,npts]  : float/double x,y,z vector

DESCRIPTION
Convert the input angular coordinate system to a 3 vector Cartesian
representation. Coordinate system axis are set up so that x points
toward the direction when priniciple angle equals 0 (ra=0,ha=0,az=0).
Y is towards increasing angle.

c1   c2    directions

ra   dec   x towards vernal equinox,y west, z celestial north.
this is a righthanded system.
ha   dec   x hour angle=0, y hour angle=pi/2(west),z=celestial north.
This is a left handed system.
az   alt   x north horizon,y east, z transit.
This is a left handed system.

vec3ToAngles
```

(See /pkg/rsi/local/libao/phil/gen/pnt/anglestovec3.pro)

## AO_AZZATORADEC_J - CONVERT AO AZ,ZA TO J2000

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
ao_azzatoradec_j - convert AO az,za to J2000
ofdate=ofdate
ARGS:
rcvnum: int    receiver number 1..16,17 for alfa,  100 for ch
use 0 to bypass model corrections.

az[n]: float/double   azimuth in degrees
za[n]: float/double za in degrees
julDat[n]: double julian date for each point.

KEYWORDS:
last[n]: double local apparent siderial time (in radians). If supplied
then this will be used rather then computing it from
the julday (you could take it from the pointing header
b.b1.h.pnt.r.lastRd
ofdate :        If set then return current ra,dec (ofdate) rather than
J2000

RETURNS:
raHrs[n]: double j2000 right ascension (hours) ..see (ofdate)
decDeg[n]: double j2000 declination (degrees)   ..see (ofdate)
DESCRIPTION:
Convert az,za encoder readings for ao to ra,dec epoch J2000. These are the
values read from the encoders after all corrections are made (model, etc).
The azimuth encoder is always the encoder on the dome side of the
azimuth arm. The routine will recompute the precession, nutation
every julian day.

EXAMPLE:
Find the ra/dec position of telescope using the az/za
positions from the .std portion of an AO header (either ri or correlator
data).

1. correlator spectral data.
suppose you have an array of correlator recs b[n]:
az=b.b1.h.std.azttd*.0001D
za=b.b1.h.std.grttd*.0001D
tm=b.b1.h.std.postmms*.001D
dayno=b.b1.h.std.date mod 1000 + tm/86400D
year =b.b1.h.std.date / 1000
2. corpwr() data.. p[n]
az=p.az*1D
za=p.za*1D
tm=p.time*1D
dayno=(p.scan/100000L) mod 1000 + tm/86400D
year =(p.scan/100000000L) + 2000L
-- warning p.scan has the last digit of the year only
after 2009 the above line doesn't work..
..The time here is the end of the integration. The az,za position
is probably from 1 sec before the end of the integration. If
you are doing 1 sec dumps then you are probably going to be
of by 1 sec of time..
3. ridata. suppose you have a long scan of ri data   d[n].
The header data is sampled once per ri record. Your accuracy will
then be the ri record spacing.
az=d.h.std.azttd*.0001D
za=d.h.std.grttd*.0001D
tm=d.h.std.postmms*.001D
dayno=d.h.std.date mod 1000 + tm/86400D
year =d.h.std.date / 1000

juldayAr=daynotojul(dayno,year) + 4.D/24.D ; +4 is ao gmt offset

The above example will have trouble with a scan that crosses midnight.
The time will go from 86399 to 0 while the day number will increment
on the next scan.
The az,za data is stored once a second. If the header records (correlator
or ri) are sampled faster than this, then you will get the value at
each second (there will be duplicates in the ra,dec array).

NOTE:
Be sure and define the juldat as a double so you have enough
precision.

```

## AO_RADECJTOAZZA - CONVERT J2000 TO AO AZ,ZA.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
ao_radecjtoazza - convert j2000 to AO az,za.
ARGS:
rcv: int    receiver number 1..12, 100 for ch
raHr[n]: float/double   J2000 ra in hours
decDeg[n]: float/double   J2000 dec in Deg
julDat[n]: double julian date for each point (utc based).
KEYWORDS:
nomodel: if set then do not include model correction in az,za.
julDat[n]: double julian date for each point (utc based).

RETURNS:
az[n]: float/double   azimuth in degrees
za[n]: float/double za in degrees
DESCRIPTION:
Convert from J2000 ra, dec to actual azimuth, zenith angle encoder
values. These are the values that would be read off of the encoder
in actual operation (the model corrections have been included). The
juliandate should be utc based. If you use julday() function remember to
added 4./24D if you use ast hour,min,seconds. The precession nutation
is computed for each julian day of the input data.

NOTE:
Be sure and define the juldat as a double so you have enough
precision.

EXAMPLE:
using the daynotojuldat routine.. you could also use juldat()
convert from hhmmss, ddmmss to hr.h deg.d
ra=205716.4D
dec=025846.1D

convert date to dayno with fraction

mon=11
day=19
year=2002
astHHMMSS=174258L
compute julianday
julday=daynotojul(daynoF,year) + 4./24.D

you could also have done this with:

julday=julday(11.D,19.D,2002.D,17.D,42.D,58.D) + 4./24.D
rcv=5               ; lbw

print,az,za
```

## AZDECTOZAHA - AZ,DECCUR TO ZA,HOUR ANGLE CONVERSION.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
azdectozaha - az,decCur to za,hour angle conversion.
SYNTAX: istat=azdectozaha(az,dec,za,ha,verbose=verbose)
ARGS:
az[n]:  double feed azimuths to use
dec[n]: double degrees declination to use. Should be current declination.

KEYWORDS:
verbose: print out the results
RETURNS:
istat  :1 if ok,-1 if at least one of the elements of az could not be computed.
eg. you picked az on wrong side of declination
za[n]  : double zenith angle in degrees
ha[n]  : double  hour angle in hours.

DESCRIPTION:
Where the equations come from:

for spherical triangle let the Angles/Sides be:
A B C the angles
a b c the opposing sides
form the triangle
a=90-dec   A=az
b=za       B=hour angle
c=90-lat   C=parallactic Angle

We are solving for b za and B hour angle.

1. The law of cosines is:
cos(a)=cos(b)cos(c) + sin(b)sin(c)cos(A)
- we need to eliminate cos(b) or sin(b) to solve for za.

2. smart page 10 formula C solve for sin(b)
sin(a)cos(C)=cos(c)sin(b) -sin(c)cos(b)cos(A)
: sin(b)=( sin(a)cos(C) + sin(c)cos(b)cos(A))/cos(c)

3. insert sin(b) into 1

cos(a)=cos(b)cos(c) + [sin(a)cos(C) + sin(c)cos(b)cos(A)]*sin(c)cos(A)/cos(c)

4. collect cos(b) terms, multiply right side by cos(c)/cos(c)

cos(b)=[cos(a)cos(c) - sin(a)cos(C)sin(c)cos(A)]
------------------------------------------
[cos^2(c) + sin^2(c)cos^2(A)]

----------------------------------------------------------------------
substitute 1-sin^2(c) for cos^2(c) -->
[1 - sin^2(c)[1-cos^2(A)]]=[1-sin^2(c)sin^2(A)]

cos(b)=[cos(a)cos(c) - sin(a)cos(C)sin(c)cos(A)]
------------------------------------------
[1-sin^2(c)sin^2(A)]
----------------------------------------------------------------------
plug in angles, arcs a,b,c from above:

cos(za)=[sin(dec)sin(lat) - cos(dec)cos(PA)cos(lat)cos(az)]
------------------------------------------
[1-cos^2(lat)sin^2(az)]

substitute PA

sin(PA)=sin(az)cos(lat)/cos(dec)   .. from law of sines
cos^2(PA)=1- sin^2(az)cos^2(lat)/cos^2(dec)
def sq:
sq = cos^2(PA)cos^2(dec)= cos^2(dec) - sin^2(az)cos^2(lat)

cos(za)=[sin(dec)sin(lat) - sqrt(sq)*cos(lat)cos(az)]
------------------------------------------
[1-cos^2(lat)sin^2(az)]
check the + and - results from sqrt(). pick za closer to 0.

this matches desh's routine that cima is using..

-->using law of sines to get hour angle
sin(ha)/sin(za)=sin(az)/sin(90-dec)
fix sign of hour angle so rising is negative

```

(See /pkg/rsi/local/libao/phil/gen/pnt/azdectozaha.pro)

## DECTOAZZA - CONVERT FROM DEC TO AZ,ZA FOR A STRIP.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
decToAzZa - convert from dec to az,za for a strip.

SYNTAX: npts=decToAzZa(dec,az,za,step=step,deg=deg,ha=ha)

ARGS:
dec : declination,degrees,minutes,seconds
unless /deg set then value is dec degrees

KEYWORDS:
step: float. seconds per step
deg :       if set then the input data in in deg, not dd mm ss
ha[n]: float if provided, then compute the az,za for the
provided hour angles (hours).&\$
zamax:   float max za to use. default:19.69

RETURNS:
npts        : int number of points returned
az[npts]    : encoder azimuth angle in degrees.
za[npts]    : encoder zenith  angle in degrees.

DESCRIPTION:
Convert from declination to azimuth, zenith angle for a
complete decstrip across the arecibo dish (latitude 18:21:14.2).
The points will be spaced step sidereal seconds in time. The routine
computes the track for +/- 2 hours and then limits the output points
to zenith angle <= 19.69 degrees.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/dectoazza.pro)

## EPVCOMPUTE - COMPUTE EARTH POS/VEL IN SOLAR SYS BARYCENTER

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
epvcompute - compute earth pos/vel in solar sys barycenter
SYNTAX: npnts=epvcompute(date,epvI,posvelI)
ARGS:
date[n]: double  mjd data and fraction of day to compute the
earth pos/velocity. The fractional part of the day
is TDB time (see times below).
epvI   : {}      info read in using epvinput. The requested dates
need to lie within the range specified in the
epvinput call.
KEYWORDS:
au     :        If set then return data as au and au/day. The
default is km and km/sec
RETURNS:
npnts:  long    the number of dates returned in posVelI.
It will return 0 if one or more of the dates requested
were not included in the epvI structure.
posVelI[6,npnts]:double the position,velocity info of the earth
center for the npnts dates requested. The order
is x,y,z,vx,vy,vz. Units are km and km/sec

DESCRPIPTION:
epvcompute returns the position and velocity of the Center of
the Earth with respect to the Solar System Barycenter at a given
TDB time.  The values are interpolated from the JPL ephemeris DE403, as
provided by their fortran routine dpleph.  Their polynomials were
re-interpolated onto approximately one-day intervals using the
numerical recipes routines.
The user inputs the daily polynomials with the epvinput() routine.
The information is stored in the epvI strucuture. The posVel info
can then be computed with this routine. The requested dates passed
into  epvcompute must lie within the range of data input by epvinput.
If not, this routine will routine 0.

The returned double posVelI[6,N] array is X, Y, Z, Vx, Vy, Vz
in J2000 equatorial coordinates (actually the IERS frame).  Units
are km and km/sec.

The format is *hard-coded* to 5th order Chebychev polynomials, and
c(1) is multiplied by 0.5, so that that needn't be done in the
evaluator.

-Mike Nolan 1997 May 29
with mods by phil 27mar98.

EXAMPLES:
;
; input the polynomial info. get the first 10 days of 2004
;
jdTomjd=2400000.5D
year=2004D
day1=1D
ndays=10
mjd1=daynotojul(day1,year) - jdtomjd
istat=epvinput(mjd1,mjd2,epvI,ndays=ndays)
;
; now compute pos/vel at 0 hours TDB start of each day.
; return data in au, au/day
;
dateAr=dblarr(10) + mjd1
n=epvcompute(dateAr,epvI,posVelI,/au)
;

NOTES:
1. The routine is vectorized so it can do multiple dates at the
same time (as long as the dates are in epvI).
2. TIMES:
TDB is barycentric dynamic time. This is the same as TDT to within
1.6 ms.
TDT tereserial dynamic time. TDT=TAI + 32.184 secs
TAI atomic time..  TAI = UTC + cumulative leap seconds (about 32
and counting).
```

(See /pkg/rsi/local/libao/phil/gen/pnt/epvcompute.pro)

## EPVINPUT - INPUT EARTH POS,VEL CHEBYCHEV POLYNOMIALS

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
epvinput - input earth pos,vel chebychev polynomials
SYNTAX: istat=epvinput(mjd1,mjd2,epvI,ndays=ndays,file=file)
ARGS:
mjd1:   long  first mjd of range to include
mdj2:   long  last mjd  of range to include (ignored if ndays is
included).
KEYWORDS:
ndays: long    if provided then ignore mjd2. The last day of range
will be mjd1+ndays-1
file : ''      If non blank then read the polynomials from here.
If null (or not provided) then use the default
filename.
RETURNS:
istat:  long    return status:
0  --> dates range is outside filename range.
> 0--> number of days returned.
< 0--> error message
-1     could not open the file.
-2     epv file header has wrong format.
-3     mjd1 is not in the file date range
epvI: {}      structure containing the polynomials. this will be
passed to epvcompute().

DESCRIPTION:
/*
epvcompute() returns the position and velocity of the Center of
the Earth with respect to the Solar System Barycenter at a given
TDB time.  The values are interpolated from the JPL ephemeris DE403, as
provided by their fortran routine dpleph.  Their polynomials were
re-interpolated onto approximately one-day intervals using the
numerical recipes routines.

epvinput() inputs the daily polynomials that are used by epvcompute()
from the file data/epv_chebfile.dat.

EPV_T1, EPV_T2 are the daily interval over which the chebychev
polynomials are interpolated.  They shouldn't necessarily be evaluated
this far out. They were chosen to be representable in base 2.
The polynomials should give correct results over the range minT,maxT.
This range is larger than 1 day. The idea is to allow some overlap in
fractional days so that you don't have to switch in the middle of an
observation.

The polynomial format is *hard-coded* to 5th order Chebychev polynomials,
and c(1) is multiplied by 0.5, so that that needn't be done in the
evaluator.

Calling sequence: call epvinput to open the data file and read the
day of interest.  Then calls to epvcompute will return the pv array.

```

(See /pkg/rsi/local/libao/phil/gen/pnt/epvinput.pro)

## GALCONV - CONVERT BETWEEN J2000 AND GALACTIC LII,BII.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
galconv - convert between J2000 and galactic LII,BII.
SYNTAX:  galconv,c1In,c2In,c1Out,c2Out,conv=conv,hms=hms,deg=deg,usened=usened
ARGS:
c1In[n] : double coord 1 In: gal l (def) or (ra) input
c2In[n] : double coord 2 In: gal b (def) or (dec) input
c1Out[n] : double coord 1 Out: ra (def) or (gal l) output
c2Out[n] : double coord 2 Out: dec (def) or (gal b) output

KEYWORDS:
togal:  By default the routine inputs galactice l,b and outputs
ra,dec. If /togal is set then the routine takes ra,dec as
input and outputs gal l,b
is set then the input parameter should be radians, and the

DESCRIPTION

Transform from  galactic l,b to J200 ra,dec. If the keyword /togal is
set then convert from ra,dec J2000 to l,b. By default the input and output
parameters are in degrees. If the /rad keyword is set then the input
and output will be in radians.

The position for the galactic pol is
NOTE:... THIS DOES NOT WORK YET. I'VE GOT 1 MORE ROTATION TO FIGURE OUT..
positions of galatic pole, center taken from ned
lb      ra           dec               ra            dec
(0,0)  266.40506655,-28.93616241 or (17:45:37.21597,-28:56:10.1847)
(0,90) 192.85949646,27.12835323  or (12:51:26.2791f, 27:07:42.0716)
```

(See /pkg/rsi/local/libao/phil/gen/pnt/galconv.pro)

## HATOAZEL - HOUR ANGLE DEC TO AZ/EL (3 VECTOR)

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
hatoazel - hour angle dec to az/el (3 vector)

ARGS:
latRd          : float/double latitude in radians

RETURNS:
azelv[3,npts]  : az,el 3 vector (for source position not feed).

DESCRIPTION
Transform from  from an hour angle dec system to an azimuth elevation system.
These are the source azimuth and elevation.

The returned 3 vector has z pointing at zenith, y pointing east (az=90),
and x pointing to north horizon (az=0). It is a left handed system.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/hatoazel.pro)

## HATORADEC - HOUR ANGLE/DEC TO RA/DEC (OF DATE).

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
hatoradec - hour angle/dec to ra/dec (of date).
ARGS:
v3[3,n]: double ha,dec data to convert to ra of date
lstRd[n]: double local apparent sidereal time in radians
OUTPUTS
v3out[3,n]:double ra,dec of date

DESCRIPTION

Transform from  from an hour angle dec system to a right ascension (of date)
dec system. The  inputs are double precision 3 vectors. The lst is passed
If the lst is the mean sidereal time , then the ra/dec and hour angle should
be the mean positions. If lst is the apparent sidereal time then the
ra/dec, hour angle should be the apparent positions.

The transformation is simlar to the raDecToHa except that we reflect
first (left handed to right handed system) before we do the rotation in
the opposite direction.
```

## JULDAYTOLMST - JULIAN DAY TO LOCAL MEAN SIDEREAL TIME.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
juldaytolmst - julian day to local mean sidereal time.
SYNTAX: lmst=juldaytolmst(juldat,obspos_deg=obspos_deg)
ARGS:
julday[n]: double julian day

KEYWORDS:
obsPos_deg: double [lat,westLong] of observatory in degrees. If
not provided then use AO position.

RETURNS :
lmst[n]: double local mean sidereal time in radians

DESCRIPTION

Convert from julian day to local mean sidereal time. By default the
latitude,long of AO is used.

If you need local apparent sidereal time, then add the equation of the
equinox to these values (see nutation_m()).

--------------------------------------------------
Mod history:
1.25feb05
there was  a bug in the code. fixed 25feb05:
An [ind] was left off of juldat0ut below:
if count gt 0 then begin
ut1Frac[ind]=ut1Frac[ind] +1.
juldat0Ut=juldat0Ut - 1L
endif
if count gt 0 then begin
ut1Frac[ind]=ut1Frac[ind] -1.
juldat0Ut=juldat0Ut + 1L
endif

- for the error to occur you had to:
1. call this routine with  juldat as an array. Calling it with
a single number would not cause the problem.
2. One of the juldats in the array had to hit 12 hrs jd
(0 hours utc).
3. the utc to ut1 correction had to be negative
Hopefully this was not a common occurence.
--------------------------------------------------
```

(See /pkg/rsi/local/libao/phil/gen/pnt/juldaytolmst.pro)

## MODEVAL - EVALUATE THE MODEL AT THE REQUESTED AZ,ZA

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
modeval - evaluate the model at the requested az,za

SYNTAX : modeval,az,za,modelData,azErrAsec,zaErrAsec,enc=enc,model=model

ARGS    :
az[]  :     azimuth positions degrees.
za[]  :     zenith angle positions degrees.
modelData: {modelData} loaded by modinp. defined in ~phil/idl/h/hdrPnt.h
azErrAsec:  [] return great circle az error in arc seconds.
zaErrAsec:  [] return great circle za error in arc seconds.

KEYWORDS:
enc   :     if 1 include encoder correction. default don't include it
model :     if 0 don't include model correction. default:include it.

DESCRIPTION:
Evaluate the model at the specified az, za locations. These are the
feed locations (not the source azimuth). Use the model data in the
structure modelData (this structure can be loaded via modinp).

Return the model errors in great circle arc seconds evaluated at the
az,za. The errors are defined such that:
1. let azComp,zaComp be the computed az, za to move the feed to.
2. compute azE, zaE from the model.
3. azTouse = azComp + AzE*asecToRad
```

(See /pkg/rsi/local/libao/phil/gen/pnt/modeval.pro)

## MODINP - INPUT MODEL COEFFICEINTS

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
modinp - input model coefficeints

SYNTAX : modinp,modelData,model=model,suf=suf,mfile=mfile,efile=efile,\$
rcv=rcv

ARGS  :
modelData: {modelData} data returned here. Structure defined in
aodefdir()/idl/h/hdrPnt.h

KEYWORDS:
model:   model name: eg modelSB, modelCB.. filename should be in
aodefdir()+"data/pnt" directory with a name modelXXX. You
can override this with the rcvNmm keyword.
default: modelSB
suf  :   suffix for the model you want. suffixes are changed as
new models are added. current suffix may00 is 11A.
if not provided then use the current model.
mfile:  string model filename (if not standard format).
efile:  string encoder filename (if not standard format)
rcvNum:  int    if supplied then use model for this receiver.
overrides model= keyword
dirmod:  string  directory for model. override default directory
WARNING:
If you are not at AO, then the model info was current when you downloaded
the aoidl distribution. It may have been updated since then. Check
the file aodefdir()+"data/pnt/lastUpdateTmStamp. It contains the date
when your data was copied from the online archive.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/modinp.pro)

## NUTATION_M - NUTATION MATRIX FOR MEAN POSITION OF DATE.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
nutation_M - nutation matrix for mean position of date.
SYNTAX: nutM=nutation_M(juldate,eqOfEq=eqOfEq)
ARGS:
juldate:    double scalar julian date for nutation.
RETURNS:
nutM[3,3]: double matrix to go from mean coordinates of date to
eqOfEq: double  if keyword supplied then the equation of the equinox will
be returned here. This is needed to go from mean to
apparent lst.
DESCRIPTION

Return nutation matrix that is used to go from mean coordinates of date to
the true coordinates of date. The nutation matrix corrects for the short
term periodic motions of the celestial pole (18 year and less....).
This matrix should be applied after the precession matrix that goes from
mean position of  epoch to mean position of date. This uses the 1980
IAU Theory of Nutation.
The equation of the equinox is also returned. This is the value to take
you from mean sidereal time to apparent sidereal time (see page B6 of AA).
It is nut (counting from 0).

Input is the julian date as a double

To create the matrix we:
1. rotate from mean equatorial system of date  to ecliptic system about x
axis (mean equinox).
2. rotate about eclipitic pole by delta psi (change in longitude due to
nutation.
3. rotate from true eclipitic system back to true equatorial by rotating
about x axis (true equinox) by -(eps + delta) eps ( mean obliquity of
eclipitic plus nutation contribution).

Since the eclipitic pole  is not affected by nutation (we ignore the
planetary contribution to the eclptic motion) only the equinox of the
eclipitic is affected by nutation. When going back from true ecliptic
to  true equatorial, use eps (avg obliquity) + deleps (obliquity due to
nutation).

Note that this method is the complete rotation versus the approximate value
used on B20 of the AE 1992.

WARNING: juldate should actually be reference to UT1 not UTC. But the
nutation matrix does not change fast enough to matter.

REFERENCE
Astronomical Almanac 1992, page B18,B20
Astronomical Almanac 1984, page S23-S26 for nutation series coefficients.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/nutation_m.pro)

## PARANGLE - COMPUTE THE PARALLACTIC ANGLE FOR AO.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
parAngle - compute the parallactic angle for AO.
SYNTAX: parAngleD=parangle(azdeg,decDeg,sitelat=sitelat)
ARGS:
azDeg  : float     the azimuth in degrees
decDeg : float     the declination of date (current) in degrees
KEYWORDS:
siteLat: float the site lattitude in degrees. If not supplied then
it uses ao's lattitude of 18:21:14.2 (but in degrees)

DESCRIPTION
Compute the parallactic angle given the source azimuth, source declination
(off date apparent) for ao. The input angles are in
degrees, the output angle is also in degrees. You can specify a different
site using the sitelatdeg

RETURNS
The parallactic angle in degrees.

REFERENCE
Spherical astronomy, smart pg 49.
Note: this probably does not yet work on arrays. need to
fix the test of  deg gt 89.9..

```

(See /pkg/rsi/local/libao/phil/gen/pnt/parangle.pro)

## PLATROTMODEL - MODEL PLATFORM ROTATION

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
platrotmodel - model platform rotation
SYNTAX: skyRot=platrotmodel(yr,lr,rotI)
ARGS:
year:   long      year for lr data (lr has a dayno but no year).
lr[n]: {}         laser ranging for dates of interest
KEYWORDS:
RETURNS:
n        : long    number of returned entries
< 0 --> error inputing model
rotI[n]: {}      array of struct holding extra tilt info.

DESCRIPTION:
A model of the platform rotations  was made for
different epochs. The models were created using data  when the platform was at the
correct height, and there was tension in all 3 tiedowns.
This routine will take laser ranging data and compute the additional
platform tilt (value - model). This could be caused by tiedowns out of position,
or because we lost tension in a tiedown.
It will also compute the pointing error caused by the excess tilt
and the pointing error caused by the current platform height.
To use this routine, you need laser ranging data with all 6 distomats
working and valid az,za and date.

The information is returned in the roti[] array of structs. It contains:
IDL> help,roti,/st
OK          LONG    1               1 rotations computed. other values are:
0 Did not fall in model data range
-1 not Enough dist measurements
-2 bad az or za value
SEC         DOUBLE  1.5369846e+09  secs from 1970 for this data point
YYMMDD      LONG    20180915       data point data as yyyymmdd (ast)
HR          FLOAT   0.169444       data point time as hr from midnight (ast)
AZ          FLOAT   95.7299        azimuth for this point (deg)
ZA          FLOAT   6.39518        za for this point (deg)
AVGH        FLOAT   382.896        avg platform height (meters). ref value is 382.896
TRINP       FLOAT   Array       plat translations dx,dy,dz in cm
ROTINP      FLOAT   Array       rotation angles [x,y] in deg for inp
ROTMOD      FLOAT   Array       rotation angles (x,y) for model at this az,za
ROTEXTR     FLOAT   Array       extra rotation (input - model) in deg about [x,y] axis.
PNTERRROT   FLOAT   0.00227990     pointing error from extra rotation (deg)
PNTERRHGHT  FLOAT   1.77322e-05    pointing error from platform height (deg)
SAVIND      LONG    2              index for model save file used for this point

NOTES:
- the x axis is toward the west
- the y axis is toward the north
- the z axis points down
- positive angles are CCW looking down the axis from the positive direction.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/platrotmodel.pro)

## PLATROTTOSKY - CONVERT PLATFORM ROTATION TO SKY ROTATION.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
platrottosky - convert platform rotation to sky rotation.
ARGS:
za[n]: float    zenith and in degrees for feed.
platRotD: float    rotation angle for platform (degrees)
KEYWORDS:
protradius: float veritcal distance from paraxial surface (za=0) where
platform rotates. default is 61 feet (close to the
connection plane of the main cables..
RETURNS:
skyRot[n]: float   Angular motion on the sky (in degrees)
retI     : {}      return a structure with the za start,za end,
inital platform angle

DESCRIPTION:

Given platform rotation angle, compute the angle moved on the sky.
The rotation angle on the sky is smaller since:
1.the radius center of curvature to feed is 435 feet.
2.the radius from center of platform rotation to feed is
60 to 100 feed (it changes as a function of za).

I've assumed the rotation is about the plane where the main cables
connect (giving a vertical distance of 61 feet paraxial surface to center
of platform rotation). I'm not sure if this is correct or not.
If you don't like that value, use the keyword protradius= to change it...

The value changes slightly with za. You can enter an array for
za and you will get back an array in skyRot.

```

(See /pkg/rsi/local/libao/phil/gen/pnt/platrottosky.pro)

## PNTERRPLATHGHT - COMPUTE POINTING ERR VS PLATFORM HEIGHT

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
pnterrplathght - compute pointing err vs platform height
SYNTAX: istat=pnterrplathght(avgHgtFt,za,pntErrAsecs,defhght=defhght,
deltahght=deltahght)
ARGS:
avgHght[n]:float  average platform heigth in feet.
za[n]: float zenith and in degrees for feed.
KEYWORDS:
defhght: float correct platform height. default is 1256.22
deltaHght:       if set then the avgHght entered is the
delta height (measured height - correct height)
RETURNS:
istat      : int   > 0 .. how many entries in pntErrAsecs
= -1   error in arguments entered
pntErrAsecs[n]:float za error in arc seconcds.

DESCRIPTION:

Compute the za pointing error when platform is moved vertically.
see http://www.naic.~phil/pnt/platMotion/pntErrVsPlatMotion.html

The default height is 1256.22 feet.
This reports the zaErr in arcseconds. The azimuth pointing error is not
affected by vertical motions.

The sign of the Error:
A height above the def hght will give a positive error. To correct for
this error you must subtract it from the model (or the requested az,za).

```

(See /pkg/rsi/local/libao/phil/gen/pnt/pnterrplathght.pro)

## PRECJ2TODATE_M - MATRIX TO PRECESS J2000 TO CURRENT.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
precj2todate_m - matrix to precess J2000 to current.
SYNTAX:  mat=precj2todate_m(juldat)
ARGS:
juldat: double   julian date in UTC system.
RETURNS:
mat[3,3]:double to precess to J2000

DESCRIPTION

Return matrix that can be used to go from J2000 equatorial rectangular
coordinates to  mean equatorial coordinates of  date (ra/dec). For
equatorial rectangular coordinates, x points to equinox, z points north,
and y points 90 to x increasing to the east.

The input consists of the UTC julian date as a double. The routine
uses the UTC time rather than converting to UT1. The difference
is not important for the usage at AO.

REFERENCE
Astronomical Almanac 1992, page B18.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/precj2todate_m.pro)

## PRECNUT - PERFORM PRECESSION AND NUTATION.

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
precNut - perform precession and nutation.
SYNTAX: precV3=precNut(vec3,juldate,toj2000=toj2000,eqOfeq=eqOfEq)
ARGS:
vec3[3,n]: double    vectors to precess
julDate  : double    julian date to use for precession matrix.

KEYWORDS:
toj2000:if set then precess from date to j2000. default is j2000 to date.

DESCRIPTION
Perform the J2000 to date or Date to J2000 precession and nutation on the
input vecotr vec3 (normalized 3 vector coordinates). The return
value is in precV3. By default the precession is from J2000 to date. The
keyword toj2000 will precess/nutate from date to J2000.

The routine computes 1 precession,nutation matrix based on the
average julDate and then applies it to all of the points n in vec3.
The juldate is normally utc based (actually it should be ut1 based but that
makes little difference here).

REFERENCE
Astronomical Almanac 1992, page B18.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/precnut.pro)

[Previous Routine] [Next Routine] [List of Routines]
```NAME:

ARGS:
raC[n] :  float current right ascension degrees
decC[n] :  float current declination degrees
lstRd[n]:  float local siderial time in radians..

KEYWORDS:

RETURNS:
az[n]    : encoder azimuth angle in degrees.
za[npts]    : encoder zenith  angle in degrees.

DESCRIPTION:
Convert from current ra,dec to az,el. this is done for the
ao dish: (latitude 18:21:14.2).
```

## RADECDIST - COMPUTE GREAT CIRCLE DISTANCE BETWEEN 2 POINTS

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
radecdist - compute great circle distance between 2 points
ARGS:
ra1[n]:  float  right ascension point 1 (hours unless /deg set)
dec1[n]:  float  declination point 1 degrees.
ra2[n]:  float  right ascension point 2 (hours unless /deg set)
dec2[n]:  float  declination point 2 degrees.
of the day.
KEYWORDS:
deg :         if set then the right ascension is in degrees, not hours

RETURNS:
ngle[n]: float  the great circle distance between the two point (in degrees).

DESCRIPTION:
Compute the great circle distance between the two ra,dec points. Return
the distance in degrees.

EXAMPLE:
ra1=11.1234     ; hours
dec1=33.345
ra2=12.321      ; hours
dec2=31.368

ra1=28.1234     ; deg
dec1=33.345
ra2=37.321      ; deg
dec2=31.368

Computation:
given points: 1,2 and let point 3 be the north pole then
the spherical astronomy cosine law is:
cosc=cosa * cosb + sina * sinb * cosC (spherical astronomy , smart pg8)

let:
A (pnt 1) angle bewteen pnt 2,3
B (pnt 2) angle between pnt 1,3
C (pnt 3) angle between 1,2 just abs(ra1-ra2)
a - distance  between b,Pole just (90- dec2)
b - distance  between a,Pole just (90 -dec1)
c - distance  between a,b .. what we want..
```

## RADECDTOHAV - CURRENT RA/DEC ANGLES TO HOURANGLE/DEC V3

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
radecDtohaV - Current ra/dec Angles to hourAngle/dec V3

ARGS  :
raRd      : float/double. right ascension of date  in radians
decRd      : float/double. declination of date in radians
lstRd[npts]: float/double. local sidereal time in radians

RETURNS:

DESCRIPTION

Transform from  from a right ascension (of date) dec system to an
hour angle dec system. The  inputs are ra,dec angles in radians and
the local sidereal time in radians.

The ra/dec system is a right handed coordinate system while the ha/dec
frame is left handed. This requires a rotation and then a reflection
to change the handedness (the minus sign around the y portion).
The returned value is the ha dec 3 vector.

RETURNS
The resulting  haDec 3 vector is returned via the pointer phaDec. The function
returns void.

```

## RADECVTOHAV - CURRENT RA/DEC (V3) TO HOURANGLE/DEC (V3)

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
radecVtohaV - current ra/dec (v3) to hourAngle/dec (v3)

ARGS  :
lstRd[npts]: float/double. local sidereal time in radians

RETURNS:

DESCRIPTION

Transform from  from a right ascension (of date) dec system to an
hour angle dec system. The  inputs are normalized 3 vectors and
that uses angles for input.

The ra/dec system is a right handed coordinate system while the ha/dec
frame is left handed. This requires a rotation and then a reflection
to change the handedness (the minus sign around the y portion).
The returned value is the ha dec 3 vector.
```

## SCHEDINFOCMP - COMPUTE SCHEDULE INFO FOR SOURCES

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
schedinfocmp - compute schedule info for sources
SYNTAX:istat=schedinfocmp(ra,dec,srcNm,srcI,coord=coord,fmt=fmt,
zamin=zamin,zamax=zamax,print=print,\$
alist=alist,hdrLabAr=hdrLabAr,exclnorise=exclnorise)
ARGS:
ra[N]: float    ra. default hms
dec[n]: ffoat   dec.. def   dms
srcNm[n]: string   source names
KEYWORDS:
coord[n]: string  coord type: j=j2000,,b=b1950,g-gal. def:j
If coord is a single element, then all elements of
ra[n],dec[n] should be this coordinate type
fmt :  int     0 - hms,dms
1= hh:mm:ss dd:mm:ss
2= degrees,degrees
3= hours,degrees
zamin : float   min za to allow
zamaz : float   max za to allow
print   :         if set the write the 1 line summaries for each line
to standard out
exclnorise:       if set then exclude sources that never rise
RETURNS:
istat  : int    n number of sources returned
-1 trouble with a source
srcI[n]: {}   struct holding info:
alist[n]: string  if supplied then return a formated list
1 line per source
hdrLabAr:string header labels for alist array.

DESCRIPTION:
Given a set of ra,decs for sources, return the  the LST rise,transit, and
set times (all in LST hours) in the structure srcI[n].
Also return the time durations: rise to transit, and keyhole to transit.
Both of these are returned in solar hours (since you'll probably want
them to compute integration times).

The routine will also return a formatted list with 1 line entry per source
if th alist=  parameter is supplied (aline[n]). the hdrLabAr is the

The user can specify the zamax,zamin where the sources rise and where they
hit the keyhole (the default is 19.69 and 1.06) .

When coordinates are precessed to the current date using the time when the
routine is run.

The structure contains:

srcI:{ srcName: ' ',\$;
raC: 0.,\$; current ra in deg
decC: 0.,\$; current dec in deg
neverRises: 0,\$; 1 (true) if never rises
;      Next 3 times are lst (sidereal) times
transitHr  : 0D,\$; transit time: lst hour
riseHr     :0D,\$ ; rises :       lst hour
setHr      :0D,\$ ; sets :        lst hour
Next 2 time durations are solar time durations
riseToTrHr :0d,\$ ; time rise to transit : solar  hours
keyHoleToTrHr:0d} ; keyhole to transit time:solar hour

```

(See /pkg/rsi/local/libao/phil/gen/pnt/schedinfocmp.pro)

## SUNAZZA - COMPUTE SUN AZ,ZA (APPROXIMATE).

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
sunazza - compute sun az,za (approximate).

SYNTAX: [az,za]=sunazza(ra1,dec1,ra2,dec2,gmtSecs)
ARGS:
ra1: double,float  hour, min,sec sun start of gmt day 1
apparent right ascension
see astronomical almanac c4->
dec1: double,float  deg, min,sec sun start of gmt day 1
ra2: double,float  hour, min,sec sun start of gmt day 2
dec2: double,float  deg, min,sec sun start of gmt day 2
gmtsecs: double        secsmidnite from ra1 day start for computation.
(ast+4 hours)
```

(See /pkg/rsi/local/libao/phil/gen/pnt/sunazza.pro)

## UTCINFOINP - INPUT PARAMETERS TO CONVERT UTC TO UT1

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
utcinfoinp - input parameters to convert utc to ut1
SYNTAX: istat=utcinfoinp(juldate,utcInfo)
ARGS:
juldate: double julian date
RETURNS:
istat  : long   0 unable to get data, 1 got the data
utcInfo:{utcInfo}

DESCRIPTION

Read the utc to ut1 information from the utcToUt1.dat file.
The UTC_INFO structure will return the information needed to go from utc to
UT1. the routine utcToUt1 converts from utc to ut1 using the information
read in here into the UTC_INFO structure. The conversion algorithm is:

utcToUt1= ( offset + ((julDay - julDayAtOffset))*rate

The offset, rate, data are input from the utcToUt1.dat file.

The user passes in the julian date and the
utcToUt1.dat file will be searched for the greatest julian date that is
less than or equal to the date passed in. If all of the values are after the
requested juliandate,  the earliest value in the file will be used and
and error will be returned.

NOTE: The file is updated whenever a leap second occurs or whenever the
drift rate changes (usually every 6 months or a year). If you have
newer versions occasionally. Check the file
aodefdir()/data/pnt/lastUpdateTmStamp for when your file was
last updated.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/utcinfoinp.pro)

## UTCTOUT1 - CONVERT UTC TO UT1

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
utctout1 - convert utc to ut1
SYNTAX: ut1FractOffset=utctout1(juldat,utcinfo)
ARGS:
julDateUtc[n]: double julian date
RETURNS:
ut1FracOffset[n]: double add this to utc based times to get ut1

DESCRIPTION

Return the offset from utc to ut1 as a fraction of a day. The returned
value (dut1Frac) is defined as ut1Frac=utcFrac + dut1Frac;
The fraction of a day can be  < 0.

The utc to ut1 conversion info is passed in via the structure UTC_INFO.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/utctout1.pro)

## UTCTOUT1NEW - CONVERT UTC TO UT1

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
utctout1new - convert utc to ut1
SYNTAX: ut1FractOffset=utctout1new(juldat)
ARGS:
julDateUtc[n]: double julian date
RETURNS:
ut1FracOffset[n]: double add this to utc based times to get ut1

DESCRIPTION

Return the offset from utc to ut1 as a fraction of a day. The returned
value (dut1Frac) is defined as ut1Frac=utcFrac + dut1Frac;
The fraction of a day can be  < 0.

The utc to ut1 conversion info is passed in via the structure UTC_INFO.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/utctout1new.pro)

## VEC3TOANGLES - CONVERT 3 VECTOR TO ANGLES

[Previous Routine] [Next Routine] [List of Routines]
```NAME:
vec3ToAngles - convert 3 vector to angles

SYNTAX: vec3toangles,v,angle1,angle2,deg=deg
ARGS:
v[3,npts]   : input 3 vector array

KEYWORDS:
c1pos       : if set then the first angle will be returned >=0
deg         : if set then return angles rather than radians

RETURNS:
c1Rd[npts]: 1st angle coordinate in radians
c2Rd[npts]: 2nd angle coordinate in radians

DESCRIPTION:
Convert the input 3 vector to angles (radians).
If posC1 is set, then the first angle will always be returned as a positive
angle (for hour angle system you may want to let the ha be negative ).

The coordinate systems are:
c1 - ra, ha, azimuth
c2 - dec,dec, altitude
```

(See /pkg/rsi/local/libao/phil/gen/pnt/vec3toangles.pro)

## VELLSRPROJ - LSR VELOCITY ALONG A GIVEN DIRECTION IN J2000 .

[Previous Routine] [List of Routines]
```NAME:
velLsrProj - lsr velocity along a given direction in J2000 .
SYNTAX: projVel=velLsrProj(dirJ2000,helioVelProj,packed=packed)
ARGS:
dirJ2000 : double 3 vector for J2000 position (unless packed keyword
set.
helioVelProj: double The observers heliocentric velocity projectd
along the dirJ2000 direction units:v/c
KEYWORDS:
packed:       if set then dirJ2000 is an array of 2 that holds
[hhmmss,ddmmss] ra and dec.
RETURNS:
projVel     : double observers lsr velocity projected onto the direction.
units: v/c. Positive velocities move away from the
coordinate center.
DESCRIPTION

Return the projected lsr velocity for an observer given the direction
in J2000 coordinates and their projected helioCentric velocity
(see velGHProj).
output:
The position for the lsr is:
ra          dec
1900  18:00:00,dec:30:00:00
J2000 18:03:50.279,30:00:16.8 from rick fischer.
raRd:4.7291355  decRd: 0.52368024

Then I converted from angles to 3vec and multiplied by 20km/sec / c

EXAMPLE:
Suppose you observered an object with 8km/sec topocentric doppler
shift, but you really wanted 8 km/sec lsr doppler shift. Here's how to
fix it (and that's why i wrote this!):
1. get the ra,dec J2000 from the pntheader.
raRd =b.b1.h.pnt.r.raJcumRd
decRd=b.b1.h.pnt.r.decJcumRd
2. convert to 3 vec
vec=anglestovec3(raRd,decRd)
3. get the observers heliocentric velocity projected onto the direction
obshelVel=b.b1.h.pnt.r.HELIOVELPROJ ; this is v/c units
4. get the projected lsr velocity of observer
obsLsrVelProj=velLsrProj(vec,obshelvel) ; this routine
5. take difference object velocity, user velocity (- since we
are interested in the relative difference of the obj,user)
vel=objVel/C - obsLsrVelProj
6. compute the doppler correction:
dopCor=(1./(1.+vel))        ; this is what should have been used.
7. Get the doppler correction used from the header
dopCorUsed=b.b1.h.dop.factor
8. compute the frequency error we made
frqErr=restFreq*(dopCor-dopCorUsed)
9. The expected line will be at frequency:dopCor*restFreq rather than
the center of the band: dopCorUsed*restFreq
cfrSbcUsed

NOTE:
To use this routine do a addpath,'gen/pnt' before calling it.
```

(See /pkg/rsi/local/libao/phil/gen/pnt/vellsrproj.pro)