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