- ABERANNUAL - COMPUTE THE ANNUAL ABERRATION OFFSET.
- AGCAZZADIF - COMPUTE AZ,ZA OFFSETS FROM AGC DATA.
- ALFABMPOS - COMPUTE ALFA BEAM POSITIONS FROM AZ,ZA,JD BEAM 0.
- ANGLESTOVEC3 - CONVERT FROM ANGLES TO 3 VECTORS.
- AO_AZZATORADEC_J - CONVERT AO AZ,ZA TO J2000
- AO_RADECJTOAZZA - CONVERT J2000 TO AO AZ,ZA.
- AZDECTOZAHA - AZ,DECCUR TO ZA,HOUR ANGLE CONVERSION.
- DECTOAZZA - CONVERT FROM DEC TO AZ,ZA FOR A STRIP.
- EPVCOMPUTE - COMPUTE EARTH POS/VEL IN SOLAR SYS BARYCENTER
- EPVINPUT - INPUT EARTH POS,VEL CHEBYCHEV POLYNOMIALS
- GALCONV - CONVERT BETWEEN J2000 AND GALACTIC LII,BII.
- HATOAZEL - HOUR ANGLE DEC TO AZ/EL (3 VECTOR)
- HATORADEC - HOUR ANGLE/DEC TO RA/DEC (OF DATE).
- JULDAYTOLMST - JULIAN DAY TO LOCAL MEAN SIDEREAL TIME.
- MODEVAL - EVALUATE THE MODEL AT THE REQUESTED AZ,ZA
- MODINP - INPUT MODEL COEFFICEINTS
- NUTATION_M - NUTATION MATRIX FOR MEAN POSITION OF DATE.
- PARANGLE - COMPUTE THE PARALLACTIC ANGLE FOR AO.
- PLATROTMODEL - MODEL PLATFORM ROTATION
- PLATROTTOSKY - CONVERT PLATFORM ROTATION TO SKY ROTATION.
- PNTERRPLATHGHT - COMPUTE POINTING ERR VS PLATFORM HEIGHT
- PRECJ2TODATE_M - MATRIX TO PRECESS J2000 TO CURRENT.
- PRECNUT - PERFORM PRECESSION AND NUTATION.
- RADECCURTOAZZA - CONVERT FROM RADECCURRENT TO AZZA
- RADECDIST - COMPUTE GREAT CIRCLE DISTANCE BETWEEN 2 POINTS
- RADECDTOHAV - CURRENT RA/DEC ANGLES TO HOURANGLE/DEC V3
- RADECVTOHAV - CURRENT RA/DEC (V3) TO HOURANGLE/DEC (V3)
- SCHEDINFOCMP - COMPUTE SCHEDULE INFO FOR SOURCES
- SUNAZZA - COMPUTE SUN AZ,ZA (APPROXIMATE).
- UTCINFOINP - INPUT PARAMETERS TO CONVERT UTC TO UT1
- UTCTOUT1 - CONVERT UTC TO UT1
- UTCTOUT1NEW - CONVERT UTC TO UT1
- VEC3TOANGLES - CONVERT 3 VECTOR TO ANGLES
- VELLSRPROJ - LSR VELOCITY ALONG A GIVEN DIRECTION IN J2000 .

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

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

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

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

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

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

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

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[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 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)**

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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