ao pointing related idl routines

Last modified: Thu Mar 16 13:07:07 2023.


List of Routines


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:
   radedJ2000 -> precNut -> raDecCurrentTrue -> addAberV->radDecCurApparent

(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.

SYNTAX: stat=agcazzadif(rcvnum,raJ,decJ,agcDat,stSec,stDay,stYear,tmStepSec,
                        npts, azdif,zadif)
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)
   read the data
       readu,lun,inparr
       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

SEE ALSO: ao_radecjtoazza

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.

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


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
added to the az,za provided.
   

(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.

SEE ALSO:
 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
SYNTAX: ao_azzatoradec_j,rcv,az,za,julDat,raHrs,decDeg,last=last
                         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
  ao_azzatoradec_j,rcvnum,az,za,juldayAr,raHrs,decDeg  

   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.

SEE ALSO:
   agcazzadif

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


AO_RADECJTOAZZA - CONVERT J2000 TO AO AZ,ZA.

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ao_radecjtoazza - convert j2000 to AO az,za.
SYNTAX: radecjtoazza,rcv,raHr,decDeg,julDat,az,za,nomodel=nomodel
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
   raH=hms1_rad(ra)*!radeg/360.D*24D
   decD=dms1_rad(dec)*!radeg

   convert date to dayno with fraction

   mon=11
   day=19
   year=2002
   astHHMMSS=174258L
   daynoF=dmtodayno(day,mon,year) + hms1_rad(asthhmss)/(2.D*!dpi)
   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

   ao_radecJtoazza,rcv,raH,decD,julday,az,za
   print,az,za

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


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,uha=uha)

ARGS:
   dec[3] : declination,degrees,minutes,seconds
            unless /deg set then value is dec[0] 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
    p12m :   int   if set then this is for the 12meter.don't switch source to 
 

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

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
     rad:  By default the input/output parameters are in degrees. If rad
           is set then the input parameter should be radians, and the
           output parameters will be radians.

 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)

SYNTAX: azelv=hatoazel(hadecV,latRd)

ARGS:
   hadecV[3,npts] : hour angle dec 3 vector (see radecdtohav).
   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).
SYNTAX: v3out=hatoradec(v3,lst)
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
 in radians.
 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.

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


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[2]: 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
    zaTouse = zaComp + ZaE*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[3] (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
                                    -3 bad date vale
 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[3]       plat translations dx,dy,dz in cm
 ROTINP      FLOAT   Array[2]       rotation angles [x,y] in deg for inp
 ROTMOD      FLOAT   Array[2]       rotation angles (x,y) for model at this az,za
 ROTEXTR     FLOAT   Array[2]       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.
SYNTAX: skyRot=platrottosky(za,platRotD,protradius=protradius,retI=retI)
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)


RADECCURTOAZZA - CONVERT FROM RADECCURRENT TO AZZA

[Previous Routine] [Next Routine] [List of Routines]
NAME:
radeccurtoazza - convert from raDecCurrent to azza

SYNTAX: radeccurtoazza,raC,decC,lstRd,az,za 

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).

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


RADECDIST - COMPUTE GREAT CIRCLE DISTANCE BETWEEN 2 POINTS

[Previous Routine] [Next Routine] [List of Routines]
NAME: 
radecdist - compute great circle distance between 2 points
SYNTAX: angle=radecdist(ra1,dec1,ra2,dec2,deg=deg)
  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
 angle=radecdist(ra1,dec1,ra2,dec2)

   ra1=28.1234     ; deg
  dec1=33.345
   ra2=37.321      ; deg
  dec2=31.368
 angle=radecdist(ra1,dec1,ra2,dec2,/deg)

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..

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


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

SYNTAX: hadecV=radecDtohaV(raRd,decRd,lstRd)

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:
   haDecV[3,npts]: hour angle/dec 3 vector.

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.

SEE ALSO: hadecVtohaV

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


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)

SYNTAX: hadecV=radecVtohaV(radecV,lstRd)

ARGS  :
 radecv[3,Npts]: float/double. 3 vector ra,dec
    lstRd[npts]: float/double. local sidereal time in radians

RETURNS:
   haDecV[3,npts]: hour angle/dec 3 vector.

DESCRIPTION

 Transform from  from a right ascension (of date) dec system to an
 hour angle dec system. The  inputs are normalized 3 vectors and
 the local sidereal time in radians. See radecDtohaV for a version
 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.

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


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[3]: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[3] is the 
header lines for these lines.

   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[3]: double,float  hour, min,sec sun start of gmt day 1
                         apparent right ascension
                    see astronomical almanac c4->
  dec1[3]: double,float  deg, min,sec sun start of gmt day 1
   ra2[3]: double,float  hour, min,sec sun start of gmt day 2
  dec2[3]: 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 
      downloaded this file from ao, then you need to redownload the
      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 
   utcInfo:{}     read in by utcinfoinp
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[3] : 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)