C routine creates a Hydra input file using a list of x and y's plus
C the header field...
	implicit real*8 (a-h,o-z)
	real*4 aweight
C data files
	character*80 filein,xyout,radecout,radectemp
	character*70 out
	parameter(maxstars=2000)
C maxstars is the maximum number of program objects that can be in the file
C it is also the maximum number of possible sky positions.
	parameter(maxguides=2000)
C "maxguides" is the maximum number of guide stars that can be in a file
	parameter(maxfibers=100)
C "maxfibers" is the maximum number of fibers in either fiber bundle
	parameter(nolimit=10000)
C "nolimit" is just a very big number
	parameter(leeslimit=4000)
C leeslimit is the maximum number of objects that can be in the list
C without Lee's software hic-uping.
C
C
C
C The different catagories are: objects, skys, and fops.
C Each has an ra, dec, and epoch
C Each has an x and y
C each has an integer id number and name
C in addition, the field center has an ra, dec, id, and name.
C (field center is at 0,0, though)
	real*8 raobjs(maxstars),decobjs(maxstars),epochobjs(maxstars),
     .         raskys(maxstars),decskys(maxstars),epochskys(maxstars),
     .      rafops(maxguides),decfops(maxguides),epochfops(maxguides)
	dimension idobjs(maxstars),idskys(maxstars),idfops(maxguides)
        character*20 nameobjs(maxstars),nameskys(maxstars),
     .               namefops(maxguides),namecent
	character*64 fieldname
	integer IDMSt(4),IDMSa(4)
	integer*4 objstatus(maxstars),skystatus(maxstars),
     .            fopstatus(maxguides),
     .            ixobjs(maxstars),iyobjs(maxstars),
     .            ixfops(maxguides),iyfops(maxguides),
     .	          ixskys(maxstars),iyskys(maxstars),
     .		  itrackfops(maxguides),
     .		  itrackskys(maxstars),itrackobjs(maxstars)
1	write(*,2)
2	format(1x,'Input file name: ',$)
	read(*,*) filein
	open(unit=1,file=filein,status='old',err=111)
	go to 3
111	write(*,6)
	go to 1
3	write(*,4)
4	format(1x,'Output root name: ',$)
	read(*,*) out
	do 8 i=1,70
	if(out(i:i).eq.' ') then
	radectemp=out(1:i-1)//'.RDTEMP'
	radecout=out(1:i-1)//'.RD'
	xyout=out(1:i-1)//'.XY'
	open(unit=8,file=radectemp,status='new',err=5)
	open(unit=9,file=radecout,status='new',err=5)
 	open(unit=10,file=xyout,status='new',err=5)
	write(*,*) xyout
	write(8,*) xyout
	go to 15
5	write(*,6)
	write(8,6)
6	format('Im sorry, but I cant open the file. Try again')
	go to 3
	end if
8	continue
	stop 'I cant handle that'
C
C
C
C Determine if this is an Old Format Nessie File
C or a New Fangled Hydra File
15	continue
	call grabvitals(fieldname,deepoch,cuepoch,aLST,exptime,
     .wave,camwave,icable,aweight,1)
C
C
C
C
C 
C
C Grab the RAs and DECs of all objects, classifying them at the
C same time
20	continue
	write(6,199)
199	format('Calling grabxy')
	call grabxy(racent,deccent,epochcent,idcent,namecent,
     .  ixfops,iyfops,epochfops,idfops,namefops,Nfops,itrackfops,
     .   maxguides,
     .  ixskys,iyskys,epochskys,idskys,nameskys,Nskys,itrackskys,
     .   maxstars,
     .  ixobjs,iyobjs,epochobjs,idobjs,nameobjs,Nobjs,itrackobjs,
     .   maxstars,
     .                 deepoch,ifile)
C bug here corrected!
	itotal=1+Nfops+Nskys+Nobjs
	if(itotal.gt.leeslimit) then
	write(6,22) itotal,leeslimit
	write(8,22) itotal,leeslimit
22	format(1x,'Total number in list=',i5, 'exceeds Hydra limit',i5)
	stop
	end if
C
	write(*,23) Nfops
	write(8,23) Nfops
23	format(1x,'There are',i5,' fops stars found')
	write(*,25) Nobjs
	write(8,25) Nobjs
25	format(1x,'There are',i5,' objects found')
C
	write(*,26) Nskys
	write(8,26) Nskys
26	format(1x,'There are',i5,' skies found')
C Now try writing them out just to see if we've got things OK:
	do 30 i=1, nobjs
	write(*,*) idobjs(i),nameobjs(i),ixobjs(i),iyobjs(i),epochobjs(i)
30	continue
	do 31 i=1, nskys
	write(*,*) idskys(i),nameskys(i),ixskys(i),iyskys(i),epochskys(i)
31	continue
	do 32 i=1, nfobs
	write(*,*) idfops(i),namefops(i),ixfops(i),iyfops(i),epochfops(i)
32	continue
	call dr2tf(3,racent,sign,IDMSt)
	call dr2af(2,deccent,sign,IDMSa)
	write(6,33) idcent,namecent, idmst,sign,idmsa,epochcent
33	format(i5,a20,1x,i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,
     &  i2.2,1x,i2.2,1x,i2.2,'.',i2.2,f7.2)
C OK now we are going to do the conversions backwards...
	angle=0.0
	write(6,123)  cuepoch,aLST,exptime,camwave,wave,angle
123	format('Calling ANTI-wiynness with',6f10.2)
	call antiwiynness(racent,deccent,epochcent,
     .               ixobjs,iyobjs,epochobjs,nobjs,raobjs,decobjs,
     .               cuepoch,aLST,exptime,camwave,wave,
     .               angle)
	write(*,123)  cuepoch,aLST,exptime,camwave,wave,angle

	call antiwiynness(racent,deccent,epochcent,
     .               ixskys,iyskys,epochskys,nskys,raskys,decskys,
     .               cuepoch,aLST,exptime,camwave,wave,
     .                angle)
	write(*,123)  cuepoch,aLST,exptime,camwave,camwave,angle

	call antiwiynness(racent,deccent,epochcent,
     .               ixfops,iyfops,epochfops,nfops,rafops,decfops,
     .               cuepoch,aLST,exptime,camwave,camwave,
     .               angle)
	call dr2tf(3,racent,sign,IDMSt)
	call dr2af(2,deccent,sign,IDMSa)
	write(*,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle
	write(8,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle
	ifred=0
	write(9,1021) namecent,IDMSt,sign,IDMSa
1021	format('   0',1x,a20,1x,
     .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,
     .' C')
Cwrite(9,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle
101	format(/1x,'Field center=',a20,
     .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,
     .i2.2,'.',i2.2,
     .2x,f7.2,' Ag=',f8.3)
	if (nobjs.lt.1) go to 110
Cwrite(9,1112)
1112	format(1x,'objects:')
	do 100 i=1,nobjs
	call dr2tf(3,raobjs(i),sign,IDMSt)
	call dr2af(2,decobjs(i),sign,IDMSa)
	write(9,1022) idobjs(i),nameobjs(i),IDMSt,sign,IDMSa

	write(8,102) idobjs(i),nameobjs(i),IDMSt,sign,IDMSa,epochobjs(i),
     .               ixobjs(i),iyobjs(i),objstatus(i)
1022	format(i4,1x,a20,1x,
     .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,
     ,' O')
102	format(i4,1x,a20,
     .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,
     ,2x,f7.2,2i7,i3)
	write(10,5000) ixobjs(i),iyobjs(i)
5000	format(2i10,' O')
100	continue
110	if(nskys.lt.1) go to 1111
Cwrite(9,103)
103	format(1x,'skys:')
	do 1101 i=1,nskys
	call dr2tf(3,raskys(i),sign,IDMSt)
	call dr2af(2,decskys(i),sign,IDMSa)
	write(8,102) idskys(i),nameskys(i),IDMSt,sign,IDMSa,epochskys(i),
     .               ixskys(i),iyskys(i),skystatus(i)
	write(9,1023) idskys(i),nameskys(i),IDMSt,sign,IDMSa
1023	format(i4,1x,a20,1x,
     .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,
     ,' S')
	write(10,5001) ixskys(i),iyskys(i)
5001	format(2i10,' S')
1101	continue
1111	if(nfops.lt.1) go to 112
	write(9,104)
104	format(1x,'fops:')
	do 105 i=1,nfops
	call dr2tf(3,rafops(i),sign,IDMSt)
	call dr2af(2,decfops(i),sign,IDMSa)
	write(8,102) idfops(i),namefops(i),IDMSt,sign,IDMSa,epochfops(i),
     .                ixfops(i),iyfops(i),fopstatus(i)
	write(10,5002) ixfops(i),iyfops(i)
	write(9,1024) idfops(i),namefops(i),IDMSt,sign,IDMSa
1024	format(i4,1x,a20,1x,
     .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,
     ,' F')

5002	format(2i10,' F')
105	continue
112	continue
	stop 'ALL DONE!'
	end
C

C
C
C
C This is a rewrite of "babyness" for WIYN
	subroutine wiynness(racent,deccent,epochcent,
     .               ra,dec,epoch,n,ix,iy,
     .		     status,
     .               cuepoch,aLST,exptime,camwave,wave,
     .               dangle)
	implicit real*8 (a-h,o-z)
	real*8 ra(*),dec(*),epoch(*)
	integer*4 ix(*),iy(*)
	integer*4 status(*)
	dimension tempmonths(12)
      	PARAMETER (AS2R=0.48481368110953599D-5) 
	data encoderscale/400./
C There are 400 encoder steps per mm of travel
	data tempmonths/33.5,33.4,35.6,41.4,49.7,57.3,61.6,60.0,57.1, 
     .                  49.9,41.0,28.3/ 
        month=(cuepoch-int(cuepoch))*12.+0.5 
        if (month.gt.12) month=month-12 
	if (month.lt.1) month=1
	if (month.gt.12) month=12
C nighttime temperature for refraction gotten from fractional part of
C current epoch
        tempadopt=tempmonths(month) 
C
C
C
C
C determine scale (encoder step per radian)
C   
C
C
C
	call mmp(racent,deccent,0D0,0D0,epochcent,cuepoch,raz,dcz)
Cwrite(9,1234) racent,deccent,epochcent
1234	format('Original  plate centers are:',2g12.6,' Epoch=',f10.4)
Cwrite(9,1233) raz,dcz,cuepoch
1233	format('Precessed plate centers are:',2g12.6,' Epoch=',f10.4)
C get current epoch ra and dec of plate center
C*****************************************
C Refraction change:
C
	iold=0
C iold=0 implies use the NEW STARLINK routines
C Get refraction coefficents A and B for current wavelen
	call getab(tempadopt,camwave,acamwave,bcamwave,iold)
	call getab(tempadopt,wave,awave,bwave,iold)
Cwrite(9,12322) tempadopt,camwave,acamwave,bcamwave
Cwrite(9,12323) tempadopt,wave,awave,bwave
12323	format('Temp=',f5.1,' At Obj wave',f7.1,' yields refract A,B coef=',
     & 2g12.6) 
12322	format('Temp=',f5.1,' At Cam wave',f7.1,' yields refract A,B coef=',
     & 2g12.6) 
C Now get refracted ra and dec of plate center
	call refractcor(raz,dcz,alst,exptime,acamwave,bcamwave,razr,
     .dczr,iold)
C razr and dczr are now the refracted, precessed centers.
C We will always use the tv wavelength for this call.
C	write(9,1235) razr,dczr
1235    format('Refracted plate centers are:',2g12.6)
C*****************************************
C
C
C then for each object:
C
C
C
C	
	do 100 i=1,n
C
C
C
C
C Do the precession thing
Cwrite(9,9100) i, ra(i),dec(i),epoch(i)
9100	format(///'Obj',i3,' Original  ra and dec=',2g12.6,
     &    ' Epoch=',f10.4)
	call mmp (ra(i),dec(i),0D0,0D0,epoch(i),cuepoch,ran,dcn)
C
C
Cwrite(9,9101) ran,dcn,cuepoch
9101	format(7x,'Precessed ra and dec=',2g12.6,' Epoch=',f10.4)
 
C
C Do the refraction thing
C***************************
C Refraction change:
C
	call refractcor(ran,dcn,alst,exptime,awave,bwave,ranr,dcnr,
     .iold)
C Note that we are using the "wave" here; may be spectrograph or
C TV depending upon subroutine call.
C
C
C Do the geometry thing on the refracted coordinates
Cwrite(9,9123) ranr,dcnr
9123	format(7x,'Refracted ra and dec=',2g12.6)
	call DS2TP(ranr,dcnr,razr,dczr,xi,eta)
Cwrite(9,123) xi,eta
 123	format('Resulting in xi, eta=',2g14.6,' (radians)')

C
C*****************************
C
C
C
C Correct for the corrector's pincushion and convert to mm.
	call correct(xi,eta,ximc,etamc)
C ximc and etamc are now in mm along the curved focal plane
Cwrite(9,9202) ximc,etamc
9202	format('Through the rotator (curved surface): ',2f10.3,' (mm)')
C
C
C
C Now determine where the projection of that would be on a plate plane.
C These are then the x and y positions before rotation and dewarping.
	call flatten(ximc,etamc,ximm,etamm)
Cwrite(9,9299) ximm,etamm
9299	format(5x,'Which projects to (flat surface):',2f10.3,' (mm)')

C
C
C Correct for rotation
	call rotate(ximm,etamm,x,y,dangle)
Cwrite(9,9203) x,y,dangle
9203	format(27x,'Rotated to:',2f10.3,' for angle=',f8.3)
C Figure out where we should plunk it down on the UNWARPED board!
	call dewarp(x,y,xp,yp)
Cwrite(9,9204) xp,yp
9204	format(22x,'Unwarped plate =',2f10.3)
C convert to encoder steps

	xpp=encoderscale*xp
	ypp=encoderscale*yp
C wow!  Now that wasn't so bad, was it!
C
C
C	write(9,9205) xpp,ypp
C9205	format(10x,'Encoder steps=',2f10.0/)
C
C
	ix(i)=xpp+0.5
	iy(i)=ypp+0.5
100	continue
	return
	end
C
C
	subroutine antiwiynness(racent,deccent,epochcent,
     .               ix,iy,epoch,n,ra,dec,
     .               cuepoch,aLST,exptime,camwave,wave,
     .               dangle)
	implicit real*8 (a-h,o-z)
	real*8 ra(*),dec(*),epoch(*)
	integer*4 ix(*),iy(*)
	dimension tempmonths(12)
      	PARAMETER (AS2R=0.48481368110953599D-5) 
	data encoderscale/400./
C There are 400 encoder steps per mm of travel
	data tempmonths/33.5,33.4,35.6,41.4,49.7,57.3,61.6,60.0,57.1, 
     .                  49.9,41.0,28.3/ 
C First, find the temperature for refraction
        month=(cuepoch-int(cuepoch))*12.+0.5 
        if (month.gt.12) month=month-12 
	if (month.lt.1) month=1
	if (month.gt.12) month=12
C nighttime temperature for refraction gotten from fractional part of
C current epoch
        tempadopt=tempmonths(month) 
C Precess plate center to current epoch:
C
C
	call mmp(racent,deccent,0D0,0D0,epochcent,cuepoch,raz,dcz)
C
C
C
        write(8,1234) racent,deccent,epochcent
1234    format('Original  plate centers are:',2g12.6,' Epoch=',f10.4)
        write(8,1233) raz,dcz,cuepoch
1233    format('Precessed plate centers are:',2g12.6,' Epoch=',f10.4)
C get current epoch ra and dec of plate center
C Refraction change:
C
	iold=0
C iold=0 implies use the NEW STARLINK routines
C Get refraction coefficents A and B for current wavelen
	call getab(tempadopt,camwave,acamwave,bcamwave,iold)
	call getab(tempadopt,wave,awave,bwave,iold)
	write(8,12322) tempadopt,camwave,acamwave,bcamwave
	write(8,12323) tempadopt,wave,awave,bwave
12323	format('Temp=',f5.1,' At Obj wave',f7.1,' yields refract A,B coef=',
     & 2g12.6) 
12322	format('Temp=',f5.1,' At Cam wave',f7.1,' yields refract A,B coef=',
     & 2g12.6) 
C Now get refracted ra and dec of plate center
	call refractcor(raz,dcz,alst,exptime,acamwave,bcamwave,razr,
     .dczr,iold)
C razr and dczr are now the refracted, precessed centers.
C We will always use the tv wavelength for this call.
	write(8,1235) razr,dczr
1235    format('Refracted plate centers are:',2g12.6)
C*****************************************
C
C
C then for each object:
C
C
C
C	
	do 100 i=1,n
C
C
C Convert to mm in floating point:
	xp=ix(i)/encoderscale
	yp=iy(i)/encoderscale
	write(8,123) ix(i),iy(i),xp,yp
123	format('ANTI-wiynness:'/'  Convesion to mm:',2i10,' --->',
     &        2f12.4)
C
C Next, go through Sam's dewarp routine...backwards
	call antidewarp(xp,yp,x,y)
	write(8,124) xp,yp,x,y
124	format('  Conversion through dewarp:',2f12.4,'-->',2f12.4)
C Now derotate these x and y's amd give me ximm and etamm

	call antirotate(x,y,ximm,etamm,dangle)
	write(8,125) x,y,ximm,etamm,dangle
125	format('Conversion through angle:',2f12.4,'--->',2f12.4,
     &   ' ang=',f10.2)
	call antiflatten(ximm,etamm,ximc,yimc)
	write(8,126) ximm, etamm,ximc,yimc
126	format('Conversion through flatten:',2f12.4,'--->',2f12.4)

C Correct for pincushion and plate scale:
	call anticorrect(ximc,yimc,xi,eta)
	write(8,127) ximc,yimc,xi,eta
127	format('Conversion through correct:',2f16.8,2g16.8)

C Go from Standard Coordinates to RA, DEC:
	call DTP2S(xi,eta,razr,dczr,ranr,dcnr)
	write(8,128) xi,eta,razr,dczr,ranr,dcnr
128	format('Converting from STD COORDS',2g16.8,' using CNT',2g16.8,
     &    '--->',2g16.8)

	call antirefract(ranr,dcnr,alst,exptime,awave,bwave,ran,dcn)
	write(8,129) ranr,dcnr,alst,exptime,awave,bwave,ran,dcn
129	format('Correcting for refraction: coords:',2g16.8/
     &  'Using LST, exptime,A, and B:',2f10.2,2g12.6/
     &  'Results in ra and dec of',2g16.8)

C Finally, let's precess;
	call mmp(ran,dcn,0D0,0D0,cuepoch,epoch(i),ra(i),dec(i))

	write(8,130) ran,dcn,cuepoch,ra(i),dec(i),epoch(i)
130	format('Precession from',2g16.8,f10.2/
     &   'Resuls in:',2g16.8)
 
100	continue
	return
	end
C
C
C
C
        subroutine antirotate(x,y,ximm,etamm,dangle)
	implicit real*8 (a-h,o-z)
C note that this angle is in degrees.
C CONVENTIONS
C remember that +eta is +N
C remember that xi is +E
C we also know that when the angle is 0, fiber 216 (0,-big) should be
C north
C so, when angle is 0, +eta (N) should be at x= 0, y=-eta*scale ----->
C  
C Using alegbra....
        etamm=x*dcosd(dangle-90.)+y*dsind(dangle-90.)
        ximm=x*dsind(dangle-90.)-y*dcosd(dangle-90.)
	return
	end
C
        subroutine rotate(ximm,etamm,x,y,dangle)
        implicit real*8 (a-h,o-z)
C note that this angle is in degrees.
C CONVENTIONS
C remember that +eta is +N
C remember that xi is +E
C we also know that when the angle is 0, fiber 216 (0,-big) should be
C north
C so, when angle is 0, +eta (N) should be at x= 0, y=-eta*scale ----->
C 
        x=etamm*dcosd(dangle-90.)+ximm*dsind(dangle-90.)
        y=-ximm*dcosd(dangle-90.)+etamm*dsind(dangle-90.)
        return
        end
C

C subroutine for getting the ra and dec out of the file
C
C
	subroutine grabxy(racent,deccent,epochcent,idcent,namecent,
     .    ixfops,iyfops,epochfops,idfops,namefops,Nfops,itrackfops,
     .     maxfops,
     .	  ixskys,iyskys,epochskys,idskys,nameskys,Nskys,itrackskys,
     .     maxskys,
     .	  ixobjs,iyobjs,epochobjs,idobjs,nameobjs,Nobjs,itrackobjs,
     .     maxobjs,
     .                 deepoch,ifile)
	implicit real*8 (a-h,o-z)
	real*8 epochobjs(*),
     .         epochskys(*),
     .         epochfops(*)
	integer*4 idobjs(*),idskys(*),idfops(*),
     .            itrackfops(*),itrackskys(*),
     .            itrackobjs(*),
     .      ixobjs(*),iyobjs(*),
     .      ixskys(*),iyskys(*),
     .      ixfops(*),iyfops(*)
        character*20 nameobjs(*),nameskys(*),
     .       namefops(*),namecent
	character*80 header
	character*1 sign
	character*1 class
C
C
C
C
C
C initialize counters
	nobjs=0
	nfops=0
	nskys=0
	icount=-1
1000	rewind(1)
C major changes made here July 18, 1991 to allow for "comment" field and
C for header field.  The following restrictions apply:
C a "header field" has a ":" somewhere in the first 17 columns
C a comment field has a # in the first column.
C everything else is assume to conform to the correct format
1001    read(1,'(A)',end=3003) header
Cwrite(*,*) header
	if(header(1:1).eq.'#') then
C comment!
        if(header(2:2).eq.'#') go to 1001
C keep the comment only if there is ONE #; ignore if there are two.
	write(9,5001) header
	write(10,5001) header
5001	format(a80)
	go to 1001
	end if
	do 33 i=1,17
	if(header(i:i).eq.':') go to 1001
33	continue
C it's a comment if there is a ":" in the first 17 characters
C
C
C OK, it's time to get serious.  First find the center.  Needs to be the first
C object.
	if(header.eq.'') go to 1001
	read(header(22:22),*,end=3003) class
Cwrite(*,*) class
	if((class.eq.'C').or.(class.eq.'c')) then
	read(header,1003) ihr,imn,sec,sign,ideg,imin,arcsec,epochcent
1003	format(34x,i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2,2x,f7.2)
	racent=STOR(ihr,imn,sec)
	deccent=ASTOR(sign,ideg,imin,arcsec)
	idcent=0
	namecent=header(15:34)
	ncent=1
	go to 1001
	end if
C Now read in each object
	read(header,1002) ix,iy,class
1002	format(2I10,1x,a1)

	icount=icount+1
C count the number of objects, staring with 0 as the first one
		if((class.eq.' ').or.(class.eq.'').or.
     &            (class.eq.'o').or.(class.eq.'O')) then
C it's an object!
		nobjs=nobjs+1
		if(nobjs.gt.maxobjs) stop 'too many objects'
		ixobjs(nobjs)=ix
		iyobjs(nobjs)=iy
		epochobjs(nobjs)=deepoch
		idobjs(nobjs)=icount
		write(nameobjs(nobjs),1004) icount
1004		format('FakeObj-',i4.4)
		itrackobjs(nobjs)=icount
	ELSE if ((class.eq.'F') .or.(class.eq.'f'))then
		nfops=nfops+1
		If (nfops.gt.maxfops) stop 'too many fops'
		ixfops(nfops)=ix
		iyfops(nfops)=iy
		epochfops(nfops)=deepoch
		idfops(nfops)=icount
		write(namefops(nfops),1005) icount
1005		format('FakeFobs-',i4.4)
		itrackfops(nfops)=icount
	ELSE if((class.eq.'S').or.(class.eq.'s')) then
		nskys=nskys+1
		If (nskys.gt.maxskys) stop 'too many skies'
		ixskys(nskys)=ix
		iyskys(nskys)=iy
		epochskys(nskys)=deepoch
		idskys(nskys)=icount
		write(nameskys(nskys),1006) icount
1006		format('FakeSkys-',i4.4)
		itrackskys(nskys)=icount
	ELSE
	
	END IF
	go to 1001
3003	Continue
	if(ncent.lt.1) stop 'cannot find the center'
	return
	end


C
C
C
C
	FUNCTION STOR(IHR,IMN,SEC)
	implicit real*8 (a-h,o-z)
	parameter(s2r=0.72722052166430398D-4)
	STOR=s2r*(60D0*(60D0*IHR+IMN)+SEC)
	return
	end
	FUNCTION ASTOR(SIGN,IDEG,IMIN,ARCSEC)
	implicit real*8 (a-h,o-z)
	character*1 SIGN
	parameter(as2r=0.48481368110953599D-5)
	if(SIGN.eq.'-') then
		ideg=-ideg
		imin=-imin
		arcsec=-arcsec
	end if
	ASTOR=as2r*(60d0*(60d0*ideg+imin)+arcsec)
	return
	end
	
		
	
	subroutine grabvitals(fieldname,deepoch,cuepoch,aLST,explength,
     .wave,camwav,icable,aweight,ifile)
C get the field name, default epoch, current epoch, LST, exp length,
C wavelength, and camera wavelength
	implicit real*8 (a-h,o-z)
	real*4 aweight
C setup the names used in the hydra file
	character*16 Things(9),ident
	character*4 color(3)
	character*80 X
	character*64 fieldname,tempname
	logical iok(7)
	character*4 cable
	character*4 weight
	character*6 wts(3)
C cheat cheat cheat cheat, to make reading in the data easy!
C in the 'C' version of course I am sure that it will be far more
C clever than this......
	data Things   /'FIELD NAME:     ', 
     .	               'INPUT EPOCH:    ', 
     .                 'CURRENT EPOCH:  ', 
     .                 'SIDEREAL TIME:  ' ,
     .                 'EXPOSURE LENGTH:' ,
     .                 'WAVELENGTH:     ',
     .                 'CABLE:          ' ,
     .                 'GUIDEWAVE:      ' ,
     .                 'WEIGHTING:      ' /
	data color/'RED ','BLUE','PINK'/
	data wts/'NONE  ','WEAK  ','STRONG'/
C global defaults:
	camwav=6000.
	icable=1
	aweight=1.0
	rewind(1)
	if(ifile.eq.0) then
	  read(1,1,end=999,err=999)
     . fieldname(1:10),deepoch,cuepoch,aLST,explength,wave
1	format(a10,5f8.2)
C format as per Nessie manual
	do 13 i=11,64
	fieldname(i:i)=' '
13	continue
	go to 100
	end if 
C Read Hydra file
11	read(1, '(a)' ,end=99) X
	read(X, '(a16)' ) ident
C If it's the old format, it's OK still.
	if(ident(1:14).eq.'SIDERIAL TIME:')
     .ident(1:14)='SIDEREAL TIME:'
	if(ident(1:11).eq.Things(1)(1:11)) then
		read(X, '(11X,a64)') fieldname
	if(fieldname(1:1).eq.' ') then
		tempname=fieldname(2:64)
		fieldname=tempname
		end if
		iok(1)=1
	end if
	if(ident(1:12).eq.Things(2)(1:12)) then
		read(X, '(12X,f12.2)') deepoch
		if((deepoch.lt.1855.).or.(deepoch.gt.2000.)) stop
     .          'trouble with input epoch'
		iok(2)=1
	end if
	if(ident(1:14).eq.Things(3)(1:14)) then
C Big Trouble if colon here!
		read(X, '(14X,f12.2)') cuepoch
		if((cuepoch.lt.1995.).or.(cuepoch.gt.2010.)) stop
     .           'trouble with current epoch'
		iok(3)=1
	end if
	if(ident(1:14).eq.Things(4)(1:14)) then
		read(X, '(14X,f12.2)') alst
		if((alst.lt.0.).or.(alst.gt.23.99999)) stop
     .               'trouble with LST entry'
		iok(4)=1
	end if
	if(ident.eq.Things(5)) then
		read(X, '(16X,f12.2)') explength
		if((explength.lt.0.).or.(explength.gt.10.)) stop
     .            'trouble with exposure length'
		iok(5)=1
	end if
	if(ident(1:11).eq.Things(6)(1:11)) then
		read(X, '(11X,f12.2)') wave
		if(wave.le.1.) wave=wave*10000.
		if((wave.gt.10000.).or.(wave.lt.3200.)) stop
     .            'trouble with wavelength'
		iok(6)=1
	end if
	if(ident(1:10).eq.Things(8)(1:10)) then
		read(X,'(10X,f12.2)') camwav
		if(camwav.le.1.) camwav=camwav*10000.
		if((camwav.gt.10000.).or.(camwav.lt.3500.)) stop
     .         'trouble with guide wavelength'
	end if
	if(ident(1:6).eq.Things(7)(1:6)) then
		read(X,'(6X,a4)') cable
		if(cable(1:1).eq.' ') then
		if((cable(2:4).ne.'RED').and.(cable(2:4).ne.'BLU').and.
     .              (cable(2:4).ne.'PIN')) stop
     .          'trouble with fiber bundle color'
		if(cable(2:4).eq.'BLU') icable=2
		if(cable(2:4).eq.'PIN') icable=3
		go to 333
		end if
	        if((cable(1:3).ne.'RED').and.(cable(1:4).ne.'BLUE').and.
     .          (cable(1:4).ne.'PINK'))stop
     .         'trouble with fiber bundle color'
		if(cable.eq.'BLUE') icable=2
		if(cable.eq.'PINK') icable=3
333		iok(7)=1
	end if
	if(ident(1:10).eq.Things(9)(1:10)) then
		read(X,'(10X,a4)') weight
		if(weight(1:1).eq. ' ') then
	    if((weight.ne.' STR').and.(weight(2:4).ne.'WEA').and.
     .      (weight(2:4).ne.'NON')) stop 'dont understand this weight'
		if(weight.eq.' STR') aweight=1.0
		if(weight(2:4).eq.'NON') aweight=0.0
		if(weight(2:4).eq.'WEA') aweight=0.5
		go to 11
		end if
		if((weight.ne.'STRO').and.(weight.ne.'WEAK').and.
     .      (weight.ne.'NONE')) stop 'dont understand this weight'
		if(weight.eq.'STRO') aweight=1.0
		if(weight.eq.'NONE') aweight=0.0
		if(weight.eq.'WEAK') aweight=0.5
	end if
		go to 11
99	do 1101 i=1,7
	if(iok(i).ne.1) then
		write(*,98) Things(i)
98		format(1x,'In new file could not find ',a16)
		stop
	end if
1101	continue
100	write(9,101) Things(1)(1:11),fieldname,
     . Things(2)(1:12),deepoch,Things(3)(1:14),cuepoch,
     . Things(4)(1:14),alst,
     . Things(5),explength,Things(6)(1:11),wave,
     . Things(7)(1:6),color(icable),Things(8)(1:10),camwav,
     . Things(9),wts(int(2.1*aweight+1))
        write(10,101) Things(1)(1:11),fieldname,
     . Things(2)(1:12),deepoch,Things(3)(1:14),cuepoch,
     . Things(4)(1:14),alst,
     . Things(5),explength,Things(6)(1:11),wave,
     . Things(7)(1:6),color(icable),Things(8)(1:10),camwav,
     . Things(9),wts(int(2.1*aweight+1))

101	format(a11,1x,a64,/a12,1x,f7.2,/a14,1x,f7.2,/a14,1x,f5.2,/a16,1x,
     .  f5.2,/
     .  a11,1x,f5.0,/a6,1x,a4,/a10,1x,f5.0,/a10,1x,a6)
	write(6,1000)
1000	format('Returning from grabvitals')
	return
999	stop 'problem reading Nessie-format input file; I quit!'
	end
C
C
C
C
C
C
C

	subroutine filetype(ifile)
C ifile=0 Nessie file
C ifile=1 Hydra file
	character*16 Thing ,ident
	data Thing  /'FIELD NAME:     '/
C Lee's file specifications claim that these two key words must be there
C in a hydra file, but that spaces are ignored and capitization is ignored;
C neither one of those things is true here, however!
C we have also cheated, in so far as we are always writing these
C keywords out as a16 and they start in column 1.
	ifile=0
	rewind(1)
C assume it is a Nessie file, unless proven otherwise!
1	read(1,2,end=99) ident
2	format(a16)
	if(ident(1:11).eq.Thing(1:11)) then
	ifile=1
	return
	end if
	go to 1
99	return
	end
C
C
C
C**********************************************
C Refraction change: refraction subroutine
	subroutine getab(temp,w,a,b,iold)
	real*8 colorfact(14)
	implicit real*8 (a-h,o-z)
	data alat/31.9633/
	parameter(as2r=0.48481368110953599D-5)
	DATA EPS/1d-10/
C precision to terminate iteration (radian)
	DATA AS/0.0065D0/
C temperature lapse rate in the troposphere (degK/meter)
	DATA UPS/.6/
C relative humidity
	DATA h0/2064./
C height of Kitt Peak
	DATA p0/795./
C typical pressure on Kitt Peak in millibars
	data colorfact/1.046,1.026,1.014,1.005,1.000,0.996,
     .0.993,0.991,0.989,0.987,0.986,0.985,0.983,0.983/
c multiply refraction at 5000A by these amounts for refraction
c at 3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,9000 A
c factor by which atmospheric pressure is reduced; i.e.,
c atmospheric pressure is 760 * P mm Hg
C check on OK wavelength w:
	if(w.lt.3000) w=3000.
	if(w.gt.10000) w=10000.
C then:
	if(iold.eq.1) then
C determine coeficients the same way as before:
	Atemp=58.3
	Btemp=-.067
C find wavelength index:		
	i=(w-3000.)/500.+1.5
	PTERM=exp(-h0/8000.)
	RTERM=colorfact(i)*PTERM/(0.962+0.0021*(TEMP-32.))
	A=Atemp*Rterm*as2r
	B=Btemp*Rterm*as2r
	else
C determine coeficients using the new & improved routines
	t0=273.+(temp-32.)*100./(212.-32.)
C turn wavelength into microns
	wl=w/10000.
C turn latitude into radians
	realat=alat*3600.*as2r
C	write(9,1235) h0,t0,p0,ups,wl,realat,as,eps
1235	format(1x,'Calling STARLINK with',8g12.6)
	call sla_REFCO(h0,t0,p0,ups,wl,realat,as,eps,a,b)
	end if
C a and b are in radians
C	write(9,1234)  w,a,b
1234	format(1x,'refraction coeficients for wavelength:',f10.2,2g12.6)
	return
	end
	subroutine antirefract(raout,decout,alst,exp,a,b,rain,decin)
C Lets just do it the silly way...
        implicit real*8 (a-h,o-z)
        inum=0
        rain=raout
        decin=decout
        write(8,*) raout,decout,alst,exp,a,b
12      call refractcor(rain,decin,alst,exp,a,b,ratest,dectest,0)
        inum=inum+1
        difra=(ratest-raout)*206265.
        difdec=(dectest-decout)*206265.
        write(8,13) inum,difra,difdec
13      format('Antirefract:',i5,2f10.2)
        if((abs(difra).lt.0.01).and.(abs(difdec).lt.0.01)) return
        rain=rain-difra/2./206265.
        decin=decin-difdec/2./206265.
        go to 12
        end

C
C
C
C
C
C This subroutine replaces "aver"
	subroutine refractcor(ra,dec,st,exptime,a,b,rar,decr,iold)
  	implicit real*8 (a-h,o-z)
	real*8 r(10),d(10)
Cwrite(9,999) ra,dec,st,exptime,a,b
999	format(1x,'Doing refractin correction for ra,dec,lst,exp,a,b=',
     ./1x,6g12.6)
	zo=dacosd(dsin(dec)*dsind(31.9633d0)+dcos(dec)*dcosd(31.9633d0)*
     .  dcos(st*15.*1.745329d-2-ra))
	division=exptime/10.
	sts=st-(4.5*division)
	SUMRA=0.
	SUMDC=0.
	SUMWEIGHTS=0.
	do 1 i=1,10
	stnow=sts+(i-1)*division
	call refract(ra,dec,stnow,a,b,rra,rdc,AIR,iold)
	r(i)=rra
	d(i)=rdc
	SUMRA=SUMRA+r(i)/AIR
	SUMDC=SUMDC+d(i)/AIR
	SUMWEIGHTS=SUMWEIGHTS+1./AIR
1	CONTINUE
	rave=sumra/sumweights
	dcave=sumdc/sumweights
C OK, now correct the RA's and Dec's
	rar=ra+rave
	decr=dec+dcave
	zp=dacosd(dsin(decr)*dsind(31.9633d0)+dcos(decr)*dcosd(31.9633d0)*
     .  dcos(st*15.*1.745329d-2-rar))
	dz=(zo-zp)*3600.
Cwrite(9,123) ra,rar,dec,decr,rave,dcave,zo,zp,dz
C23	format(1x,'REFRACTION ra,dec=',
C    .1x,4g12.6,/1x,2g12.4,2f10.4,f6.1)
	return
	end
      subroutine refract(ra,dec,st,a,b,rra,rdc,air,iold)
c INPUT:
c ra is ra in radians
c dec is dec in radians
c st is siderial time in decimal hrs
C a and b are the refraction coefficients for r=A*tanz + Btan^3z
c OUTPUT:
c rra is correction for refraction in right ascention in RADIANS
c                    corrected for declination
c rdc is correction for refraction in declination in RADIANS
c air is the airmass
	implicit real*8 (a-h,o-z)
	DATA degtorad/0.0174532925d0/
	data pi/3.141592653589793d0/
	data alat/31.9633/
	realat=alat*pi/180.
	twopi=pi*2.
	RST=ST*15.*pi/180.
	HAold=RST-RA
	if(haold.ge.pi)haold=haold-twopi
	if(haold.lt.-pi)haold=haold+twopi
C see to it that ha is between -pi and +pi
	cz=dsin(dec)*dsin(realat)+dcos(dec)*dcos(realat)*dcos(HAold)
	air=1./cz
	z=dacos(cz)
Cwrite(9,991) ra,dec,st,a,b,z,air
C91	format(1x,'unrefracted ra,dec,st,a,b,z,air=',
C    .6f12.4,f10.2)
C zold is now the old zenith distance
	if(iold.eq.1) then
	zr=z-A*dtan(z)+B*(dtan(z))**3.
C this is almost right.  A better thing would be to do
	dz=(z-zr)*206280.d0
Cwrite(9,7191) z,zr,dz
C191	format(1x,'With OLDCODE: z, zr, dz=',2g12.6,f10.2)
	else
	call sla_refz(z,a,b,zr)
	dz=(z-zr)*206280.d0
Cwrite(9,7192) z,zr,dz
C192	format(1x,'With STARLINK: z, zr, dz=',2g12.6,f10.2)
	end if
C znew is now the new zenith distance
Cwrite(9,990) zr
C90	format(1x,'Leading to new zenith distance=',f12.4)
	x=dsin(HAold)/(dcos(HAold)*dsin(realat)-dtan(dec)*dcos(realat))
C x is now the arctan of the azimuth:
	if(Haold.ge.0.) then
		if(x.ge.0.0d0) then
	        Az=pi+datan(x)
		else
		Az=twopi+datan(x)
		end if
	else
		if(x.ge.0.d0) then
		Az=datan(x)
		else
		Az=pi+datan(x)
		end if
	end if
811	format(1x,'z,zr,HAold,decold,decnew,x,Az=',
     ./7g12.6)
C Az should be between 0 and twopi at this point.
C find new declination:
	s=dsin(realat)*dcos(zr)+dcos(realat)*dsin(zr)*dcos(Az)
	decnew=dasin(s)
Cwrite(9,811)z,zr,Haold,dec,decnew,x,Az
	rdc=decnew-dec
C find new hour hangle:
	sinha=-dsin(Az)*dsin(zr)/dcos(decnew)
	cosha=(dcos(zr)-dsin(decnew)*dsin(realat))/
     .                  (dcos(decnew)*dcos(realat))
	hanew=datan2(sinha,cosha)
	rra=haold-hanew
	if(rra.gt.pi) rra=rra-twopi
	if(rra.lt.-pi) rra=rra+twopi
Cwrite(9,780) rra,rdc
C80	format(1x,'Refraction=',2g12.6)
C change has to be between -pi and pi (and of course is usually much smaller!)
	return
	end
        subroutine anticorrect(ximc,etamc,xi,eta)
        implicit real*8 (a-h,o-z)
*  ARC SECONDS TO RADIANS
         PARAMETER (AS2R=4.8481368110953599D-6)
Cdata scale/9.4032/
        data scale/9.4000/
C xout is in mm
        xout=ximc
        yout=etamc
C xi and eta are in radians
        xi=xout*scale*as2r
        eta=yout*scale*as2r
C Start with the assumption that the output is similar to the input
C xi and eta are in radians; xtest and ytest are in mm
12      call correct(xi,eta,xtest,ytest)
        write(6,*) xi,eta,xtest,ytest
        inum=inum+1
C difx is in mm
        difx=xtest-xout
        dify=ytest-yout
        if ((abs(difx).lt. 0.00001).and.(abs(dify).lt. 0.00001)) return
C Take half the difference
C Put it back in radians:
        xi=xi-(difx/2.)*scale*as2r
        eta=eta-(dify/2.)*scale*as2r
        go to 12
        end

	subroutine correct(xi,eta,ximc,etamc)
C this program considerably simplied after our initial astrometry 3/95
C old code:
C	implicit real*8 (a-h,o-z)
C	data AA,BB,CC,DD/6.391006d0,1.249909d-4,2.576783d-5,4.136236d-7/
C	d=dsqrt(xi**2+eta**2)
C	dp=AA*d+BB*d**2+CC*d**3+DD*d**4
C Great.  Now, what are the new xi's and eta's?
C        ximc=xi*(dp/d)
C	 etamc=eta*(dp/d)
C ximc and etamc are now xi and eta in the unlikely units of "mm"
C
	implicit real*8 (a-h,o-z)
*  ARC SECONDS TO RADIANS   
      PARAMETER (AS2R=4.8481368110953599D-6)

Cdata scale/9.4032/ 
	data scale/9.4000/
C this done for Sam 8/28/95
	data c/81.0/
C Correct for pincushion distortion:
	f=1.+c*(xi**2+eta**2)
	xip=xi*f
	etap=eta*f
Cwrite(9,101) f,xip,etap
101	format('Pincushion factor=',f10.6/
     & 'Resulting in xi, eta=',2g14.6)
C Correct for plate scale
	ximc=xip/scale/as2r
	etamc=etap/scale/as2r
C OK: ximc and etamc are now eta and xi in mm along curved surface
	return
	end
C
C
	subroutine antidewarp(xout,yout,xin,yin)
	implicit real*8 (a-h,o-z)
C Cheat!  Just keep calling dewarp until we get good answers:
	inum=0
	xin=xout
	yin=yout
C Start with the assumption that xin and yin are pretty close to
C to the final answers:
12	call dewarp(xin,yin,xtest,ytest)
	inum=inum+1
	difx=xtest-xout
	dify=ytest-yout
	write(8,13) inum,difx,dify
13	format('antidewarp:',i5,2f10.3)
	if ((abs(difx).lt.0.0001).and.(abs(dify).lt.0.0001)) return
C Take half the difference
	xin=xin-difx/2.
	yin=yin-dify/2.
	go to 12
	end

C warp subroutine:
	subroutine dewarp(xin,yin,xout,yout)
	implicit real*8 (a-h,o-z)
C These are Sam's estimage of the radius from the pupil to the focal
C surface
	data rfocal/5675.2/

C
C What is the actual orthogonal distance that we want from the
C radius vector to the top of the fiber.
	din=dsqrt(xin**2+yin**2)
C The following were modified as per Sam's 13 March note
C re the "mid radius" where length is conserved....

	if(din.lt.151.9325) then
	reff=5679.881
	else
	temp=din*400.
	reff=5670.799 + temp*(2.988726e-4 - 2.458933e-9 * temp)
C Note that Sam's equation was in terms of encoder units.
	end if
C Great.  Now, what is the distance we want along the plate prior to
C being warped?
	dout=reff*dasin(din/rfocal)
C Now make sure that the proportions work out:
	xout=0.
	yout=0.
	if(din.gt.0.) then
	yout=(dout/din)*yin
	xout=(dout/din)*xin
	end if
Cwrite(9,1234) din,dout,reff,rfocal,xin,xout,yin,yout
1234	format('Din='f10.3,' Dout=',f10.3,' reff=',f10.3, 'rfoc=',
     &  f10.3,'xin,out; yin,out=',f10.3)
	return
	end
C Anti-deflatten!
	        subroutine antiflatten(xout,yout,xin,yin)
        implicit real*8 (a-h,o-z)
        xin=xout
        yin=yout
C Start with the assumption that the output is similar to the input
12      call flatten(xin,yin,xtest,ytest)
        inum=inum+1
        difx=xtest-xout
        dify=ytest-yout
        write(8,13) inum,difx,dify
13      format('antiflatten:',i5,2f10.3)
        if ((abs(difx).lt.0.0001).and.(abs(dify).lt.0.0001)) return
C Take half the difference
        xin=xin-difx/2.
        yin=yin-dify/2.
        go to 12
        end
C
C This subroutine re-written and considerably simplied from old Hydra
C days!
	subroutine flatten(x,y,xp,yp)
	implicit real*8 (a-h,o-z)
	r=5675.2
C = 223.4 inches radius of curvature
c r is radius of curve (mm) of Lockheed camera plate adapted for WIYN  
	zed=sqrt(x**2+y**2)
	theta=zed/r
	h=r*sin(theta)
	xp=0.
	yp=0.
	if (zed.gt.0.) then
	factor=h/zed
	xp=x*factor
	yp=y*factor
	end if
	return
	end
C
C subroutines used (from kaboomlib.a)
C
C call libinit() [no arguments]: initializes the "data base"; use before any
C other of the C calls
C 
C call rotmode (flag) turns rotation on (1) or off (0).
C
C call rotangle (real*8 angle in radians): sets the rotation angle
C 
C call concenmode (flag) turns concentricity correction on (1) or off (0)
C
C call setfib (ibut, x,y,istat) sets button ibut to integer "sky" coordinates
C x,y.  Returns istat (integer) 0 = sucessful; 1 = if fiber would be over-
C extended or pushed back beyond the park circle; 2= if the attempted move was
C to a step
C
C call checkfib(ibut,x,y,istat) checks to see if button ibut can go to position
C x,y; same as setfib but the database is not updated.
C 
C call setfibphys/checkfibphys [same arguments as above]: does the same but
C uses PHYSICAL rather than SKY coordinates.  I use this if whereis has
C returns the x and y position of a button that I may need to put back
C there.
C 
C call whereis (ibut,x1,y1,x2,y2,x3,y3) : returns the location (x1,y1) of
C button ibut [ PHYSCIAL coordinates].  x2,y2 should always be the same
C as (x1,y1), and are the coordinates used for collision testing
C (x3,y3) are the taget location for the next move and are only sometimes
C valid
C 
C call park (ibut,flag) Parks the fiber; 0 = park; 1 = stow position
C
C ouch=interract (ibuta,ibutb) 0 = they aren't in each other's personal
C space; 1=they are
C 
C ouch=ibuttonhitsanother returns 0 if there is no interraction with
C another button; !=0 if there is.  I don't use this.
C
C the rest of these stuff is directly from the STARLINK library
      SUBROUTINE sla_REFRO (Z0,H0,T0,P0,UPS,WL,PH,AS,EPS,REF)

*+
*
*     - - - - - -
*      R E F R O
*     - - - - - -
*
*
*  Atmospheric refraction for optical wavelengths
*
*
*  Given:
*    Z0    dp    observed zenith distance of the source (radian)
*    H0    dp    height of the observer above sea level (metre)
*    T0    dp    ambient temperature at the observer (deg K)
*    P0    dp    pressure at the observer (millibar)
*    UPS   dp    relative humidity at the observer (range 0-1)
*    WL    dp    effective wavelength of the source (micrometre)
*    PH    dp    latitude of the observer (radian, astronomical)
*    AS    dp    temperature lapse rate in the troposphere (degK/metre)
*    EPS   dp    precision required to terminate iteration (radian)
*
*  Returned:
*    REF   dp    refraction: in vacuo ZD minus observed ZD (radian)
*
*  Called:  sla_$ATMT, sla_$ATMS
*
*  Typical values for the AS and EPS arguments would be 0.0065D0 and
*  1D-10 respectively.
*
*  This routine computes the refraction for zenith distances up to
*  and beyond 90 deg using the method of Hohenkerk and Sinclair (NAO
*  Technical Notes 59 and 63).  The code is a trivially modified form
*  of the AREF subroutine of C.Hohenkerk (HMNAO, September 1984).
*  Most of the modifications are cosmetic;  in addition the angle
*  arguments have been changed to radians and the Z0 = 0 case is
*  handled.  For Z0 > 90 deg it is up to the caller to supply an
*  H0 consistent with the depression of the horizon thus implied.
*
*  Fixed values of the water vapour exponent, height of tropopause, and
*  height at which refraction is negligible are used.
*
*  Original version by C.Hohenkerk, HMNAO, September 1986.
*  Starlink version by P.T.Wallace, July 1986.
*
*+

      IMPLICIT NONE

      DOUBLE PRECISION Z0,H0,T0,P0,UPS,WL,PH,AS,EPS,REF

*
*  Fixed parameters
*
      DOUBLE PRECISION GCR,DMD,DMW,S,GAMMA,HT,HS
*  Universal gas constant
      PARAMETER (GCR=8314.36D0)
*  Molecular weight of dry air
      PARAMETER (DMD=28.966D0)
*  Molecular weight of water vapour
      PARAMETER (DMW=18.016D0)
*  Mean Earth radius (metre)
      PARAMETER (S=6378120D0)
*  Exponent of temperature dependence of water vapour pressure
      PARAMETER (GAMMA=18.36D0)
*  Height of tropopause (metre)
      PARAMETER (HT=11000D0)
*  Upper limit for refractive effects (metre)
      PARAMETER (HS=80000D0)

      INTEGER IN,IS,K,ISTART,I,J
      DOUBLE PRECISION WLSQ,GB,Z1,A(10),PW0,R0,T0O,DN0,DNDR0,
     :                 SK0,F0,RT,TT,DNT,DNDRT,ZT,FT,DNTS,DNDRTS,
     :                 SINE,ZTS,FTS,RS,DNS,DNDRS,ZS,FS,REFO,
     :                 FE,FO,H,FB,FF,STEP,Z,R,RG,TG,
     :                 T,DN,DNDR,F,REFP,REFT

*  The refraction integrand
      DOUBLE PRECISION REFI
      REFI(R,DN,DNDR) = R*DNDR/(DN+R*DNDR)



*  Set up model atmosphere parameters defined at the observer
      WLSQ = WL*WL
      GB= 9.784D0*(1D0-0.0026D0*COS(2D0*PH)-0.00000028D0*H0)
      Z1 = ((287.604D0+1.6288D0/WLSQ+0.0136D0/(WLSQ*WLSQ))
     :         *273.15D0/1013.25D0)*1D-6
      A(1) = ABS(AS)
      A(2) = (GB*DMD)/GCR
      A(3) = A(2)/A(1)
      A(4) = GAMMA
      PW0 = UPS*(T0/247.1D0)**A(4)
      A(5) = PW0*(1D0-DMW/DMD)*A(3)/(A(4)-A(3))
      A(6) = P0+A(5)
      A(7) = Z1*A(6)/T0
      A(8) = (Z1*A(5)+11.2684D-6*PW0)/T0
      A(9) = (A(3)-1D0)*A(1)*A(7)/T0
      A(10) = (A(4)-1D0)*A(1)*A(8)/T0
      A(3) = A(3)-2D0
      A(4) = A(4)-2D0

*  At the observer
      R0 = S+H0
      CALL sla_$ATMT(R0,T0,A,R0,T0O,DN0,DNDR0)
      SK0 = DN0*R0*SIN(Z0)
      F0 = REFI(R0,DN0,DNDR0)

*  At the tropopause in the troposphere
      RT = S+HT
      CALL sla_$ATMT(R0,T0,A,RT,TT,DNT,DNDRT)
      SINE = SK0/(RT*DNT)
      ZT = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0)))
      FT = REFI(RT,DNT,DNDRT)

*  At the tropopause in the stratosphere
      CALL sla_$ATMS(RT,TT,DNT,A(2),RT,DNTS,DNDRTS)
      SINE = SK0/(RT*DNTS)
      ZTS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0)))
      FTS = REFI(RT,DNTS,DNDRTS)

*  At the stratosphere limit
      RS = S+HS
      CALL sla_$ATMS(RT,TT,DNT,A(2),RS,DNS,DNDRS)
      SINE = SK0/(RS*DNS)
      ZS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0)))
      FS = REFI(RS,DNS,DNDRS)

*
*  Integrate the refraction integral in two parts;  first in the
*  troposphere (K=1), then in the stratosphere (K=2)
*
      REFO = -999.999D0
      IS = 16
      DO K = 1,2
         ISTART = 0
         FE = 0D0
         FO = 0D0

         IF (K.EQ.1) THEN
            H = (ZT-Z0)/DBLE(IS)
            FB = F0
            FF = FT
         ELSE
            H = (ZS-ZTS)/DBLE(IS)
            FB = FTS
            FF = FS
         ENDIF

         IN = IS-1
         IS = IS/2
         STEP = H
*
*     Start of iteration loop (terminates when specified precision reached)
 200     CONTINUE

         DO I=1,IN

            IF (I.EQ.1.AND.K.EQ.1) THEN
               Z = Z0 + H
               R = R0
            ELSE IF (I.EQ.1.AND.K.EQ.2) THEN
               Z = ZTS + H
               R = RT
            ELSE
               Z = Z + STEP
            ENDIF

*        Given the zenith distance (Z) find R
            RG = R
            DO J=1,4
               IF (K.EQ.1) THEN
                  CALL sla_$ATMT(R0,T0,A,RG,TG,DN,DNDR)
               ELSE
                  CALL sla_$ATMS(RT,TT,DNT,A(2),RG,DN,DNDR)
               ENDIF
               IF (Z.GT.1D-20)
     :               RG = RG-((RG*DN-SK0/SIN(Z))/(DN+RG*DNDR))
            ENDDO
            R = RG

*        Find refractive index and integrand at R
            IF (K.EQ.1) THEN
               CALL sla_$ATMT(R0,T0,A,R,T,DN,DNDR)
            ELSE
               CALL sla_$ATMS(RT,TT,DNT,A(2),R,DN,DNDR)
            ENDIF

            F = REFI(R,DN,DNDR)

            IF (ISTART.EQ.0.AND.MOD(I,2).EQ.0) THEN
               FE = FE+F
            ELSE
               FO = FO+F
            ENDIF
         ENDDO

*     Evaluate the integrand using Simpson's Rule
         REFP = H*(FB+4D0*FO+2D0*FE+FF)/3D0

*     Skip if required precision reached
         IF (ABS(REFP-REFO).LE.0.5D0*EPS) GO TO 900

*     Prepare for next iteration
         IS = 2*IS
         IN = IS
         STEP = H
         H = H/2D0
         FE = FE+FO
         FO = 0D0
         REFO = REFP
         IF (ISTART.EQ.0) ISTART=1

*     Iterate
         GO TO 200

*     Save troposphere component
 900     CONTINUE
         IF (K.EQ.1) REFT=REFP 

      ENDDO

*  Result
      REF=REFT+REFP

      END
      SUBROUTINE sla_$ATMT (R0,T0,A,R,T,DN,DNDR)

*+
*
*     - - - - -
*      A T M T
*     - - - - -
*
*
*  Internal routine used by REFRO
*
*  Model atmosphere for the troposphere
*
*
*  Given:
*    RT    dp    height of tropopause from centre of the Earth (metre)
*    R0    dp    height of observer from centre of the Earth (metre)
*    T0    dp    temperature at the observer (deg K)
*    A   dp(10)  constants defined at the observer
*    R     dp    current distance from the centre of the Earth (metre)
*
*  Returned:
*    T     dp    temperature at R (deg K)
*    DN    dp    refractive index at R
*    DNDR  dp    rate the refractive index is changing at R
*
*  This routine is a version of the ATMOSTRO routine (C.Hohenkerk,
*  HMNAO), with trivial modifications.
*
*  Original version by C.Hohenkerk, HMNAO, August 1984
*  Starlink version by P.T.Wallace, July 1986
*
*+

      IMPLICIT NONE

      DOUBLE PRECISION R0,T0,A(10),R,T,DN,DNDR

      DOUBLE PRECISION TT0,TT01,TT02


      T = T0-A(1)*(R-R0)
      TT0 = T/T0
      TT01 = TT0**(A(3))
      TT02 = TT0**(A(4))
      DN = 1D0+(A(7)*TT01-A(8)*TT02)*TT0
      DNDR = -A(9)*TT01+A(10)*TT02

      END
      SUBROUTINE sla_$ATMS (RT,TT,DNT,GMDR,R,DN,DNDR)

*+
*
*     - - - - -
*      A T M S
*     - - - - -
*
*
*  Internal routine used by REFRO
*
*  Model atmosphere for the stratosphere
*
*
*  Given:
*    RT    dp    height of tropopause from centre of the Earth (metre)
*    TT    dp    temperature at the tropopause (deg K)
*    DNT   dp    refractive index at the tropopause
*    GMDR  dp    constant of the atmospheric model = G*MD/R
*    R     dp    current distance from the centre of the Earth (metre)
*
*  Returned:
*    DN    dp    refractive index at R
*    DNDR  dp    rate the refractive index is changing at R
*
*  This routine is a version of the ATMOSSTR routine (C.Hohenkerk,
*  HMNAO), with trivial modifications.
*
*  Original version by C.Hohenkerk, HMNAO, August 1984
*  Starlink version by P.T.Wallace, July 1986
*
*+

      IMPLICIT NONE

      DOUBLE PRECISION RT,TT,DNT,GMDR,R,DN,DNDR

      DOUBLE PRECISION B



      B = GMDR/TT
      DN = 1D0+(DNT-1D0)*EXP(-B*(R-RT))
      DNDR = -B*(DNT-1D0)*EXP(-B*(R-RT))

      END
      SUBROUTINE sla_REFZ (ZU,REFA,REFB,ZR)

*+
*
*     - - - - -
*      R E F Z
*     - - - - -
*
*
*  Adjust an unrefracted zenith distance to include the effect of
*  atmospheric refraction, using the simple A tan Z + B tan**3 Z
*  model.
*
*
*  Given:
*    ZU    dp    unrefracted zenith distance of the source (radian)
*    REFA  dp    tan Z coefficient (radian)
*    REFB  dp    tan**3 Z coefficient (radian)
*
*  Returned:
*    ZR    dp    refracted zenith distance (radian)
*
*  Note that this routine applies the adjustment for refraction in
*  the opposite sense to the usual one - it takes an unrefracted
*  (in vacuo) position and produces an observed (refracted)
*  position, whereas the basic A tan Z + B tan**3 Z model strictly
*  applies to the case where a refracted position is available and
*  must be corrected for refraction.  This requires an inverted form of
*  the refraction expression; the formula used here is taken from
*  RGO/LaPalma Technical Note No. 42, "A control system for
*  the William Herschel Telescope", by R.Laing.
*
*  For ZD > 85 deg, a fixed value ZD = 85 deg is used.
*
*  See also the routine sla_REFV, which performs the adjustment in
*  Cartesian Az/El coordinates.
*
*  P.T.Wallace   Starlink   July 1986
*
*+

      IMPLICIT NONE

      DOUBLE PRECISION ZU,REFA,REFB,ZR

      DOUBLE PRECISION ZL,S,C,T,TSQ,TCU



*  Keep ZD within usable range
      ZL = MIN(ABS(ZU),1.483529864D0)

*  Functions of ZD
      S = SIN(ZL)
      C = COS(ZL)
      T = S/C
      TSQ = T*T
      TCU = T*TSQ

*  Refracted ZD
      ZR = ZU-(REFA*T+REFB*TCU)/(1D0+(REFA+3D0*REFB*TSQ)/(C*C))

      END
      SUBROUTINE sla_REFCO (H0,T0,P0,UPS,WL,PH,AS,EPS,REFA,REFB)

*+
*
*     - - - - - -
*      R E F C O
*     - - - - - -
*
*
*  Determine constants A and B in atmospheric refraction model
*  dZ = A tan Z + B tan**3 Z (for optical wavelengths).
*
*
*  Given:
*    H0    dp    height of the observer above sea level (metre)
*    T0    dp    ambient temperature at the observer (deg K)
*    P0    dp    pressure at the observer (millibar)
*    UPS   dp    relative humidity at the observer (range 0-1)
*    WL    dp    effective wavelength of the source (micrometre)
*    PH    dp    latitude of the observer (radian, astronomical)
*    AS    dp    temperature lapse rate in the troposphere (degK/metre)
*    EPS   dp    precision required to terminate iteration (radian)
*
*  Returned:
*    REFA  dp    tan Z coefficient (radian)
*    REFB  dp    tan**3 coefficient (radian)
*
*  Called:  sla_REFRO
*
*  Typical values for the AS and EPS arguments would be 0.0065D0 and
*  1D-10 respectively.
*
*  Relative to the comprehensive refraction model used by this routine,
*  the simple A tan Z + B tan**3 Z formula achieves 0.5 arcsec accuracy
*  for ZD < 80 deg, 0.01 arcsec accuracy for ZD < 60 deg, and
*  0.001 arcsec accuracy for ZD < 45 deg.
*
*  P.T.Wallace   Starlink   July 1987
*
*+

      IMPLICIT NONE

      DOUBLE PRECISION H0,T0,P0,UPS,WL,PH,AS,EPS,REFA,REFB

      DOUBLE PRECISION ATN1,ATN4,R1,R2

*  Sample zenith distances: arctan(1) and arctan(4)
      PARAMETER (ATN1=0.7853981634D0,
     :           ATN4=1.325817664D0)



*  Determine refraction for the two sample zenith distances
      CALL sla_REFRO(ATN1,H0,T0,P0,UPS,WL,PH,AS,EPS,R1)
      CALL sla_REFRO(ATN4,H0,T0,P0,UPS,WL,PH,AS,EPS,R2)

*  Solve for refraction constants
      REFA = (64D0*R1-R2)/60D0
      REFB = (R2-4D0*R1)/60D0

      END
        SUBROUTINE RDCAL(X,Y,PCRA,PCDEC,ITEL,RA,DEC)  
    
*  ADJUST IDEAL PLATE COORDINATES (X,Y) ACCORDING TO
*  TELESCOPE TYPE GIVING TANGENTIAL COORDS (XI,ETA) 
*  THEN CALCULATES (RA,DEC) 
    
    
*  GIVEN:   
*     X       IDEAL PLATE COORDINATE
*     Y       IDEAL PLATE COORDINATE
*     PCRA    PLATE CENTRE RA   
*     PCDEC   PLATE CENTRE DEC  
*     ITEL    TELESCOPE TYPE
    
*  RETURNED:
*     RA      CALCULATED RA 
*     DEC     CALCULATED DEC
    
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION X,Y,PCRA,PCDEC   
      INTEGER ITEL  
      DOUBLE PRECISION RA,DEC   
    
    
*  ADJUST FOR TELESCOPE TYPE
      GO TO (11,12,13,14),ITEL  
    
*  ASTROGRAPH   
   11 CONTINUE  
      XI=X  
      ETA=Y 
      GO TO 20  
    
*  SCHMIDT  
   12 CONTINUE  
      CALL SCHMIT(X,Y,CURV) 
      XI=X*CURV 
      ETA=Y*CURV
      GO TO 20  
    
*  AAT DOUBLET  
   13 CONTINUE  
      XI=X  
      ETA=Y 
      CALL PNCUSH(ITEL,XI,ETA,'DU') 
      GO TO 20  
    
*  AAT TRIPLET  
   14 CONTINUE  
      XI=X  
      ETA=Y 
      CALL PNCUSH(ITEL,XI,ETA,'DU') 
   20 CONTINUE  
    
*  COMPUTE (RA,DEC) 
      CALL DTP2S(XI,ETA,PCRA,PCDEC,RA,DEC)  
    
      END   
      SUBROUTINE PNCUSH(ITEL,X,Y,KDIRN) 
    
*  ALLOW FOR PINCUSHION DISTORTION  
    
*  GIVEN:   
*     ITEL    TELESCOPE TYPE : 3 = AAT2, 4 = AAT3   
*     X,Y     PLATE COORDINATES (RADIANS)   
*     KDIRN   DIRECTION OF ADJUSTMENT   
    
*  RETURNED:
*     X,Y     ADJUSTED  
    
    
*  KDIRN VALUES:
*     'UD'    FROM TANGENT-PLANE COORDINATES
*             PREDICTS DISTORTED COORDINATES
*     'DU'    FROM DISTORTED COORDINATES
*             PREDICTS TANGENT-PLANE COORDINATES
    
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      INTEGER ITEL  
      DOUBLE PRECISION X,Y  
      CHARACTER*2 KDIRN 
    
*  PINCUSHION DISTORTION COEFFICIENTS : AAT DOUBLET, TRIPLET
*     (DOUBLET USED AS ERROR DEFAULT)   
    
      PARAMETER (PCAAT2=147.069D0)  
      PARAMETER (PCAAT3=178.585D0)  
    
    
*  PICK UP CORRECT DISTORTION COEFFICIENT   
      D=PCAAT3  
      IF(ITEL.NE.4) D=PCAAT2
    
*  DISTANCE FROM CENTRE 
      RSQ=X*X + Y*Y 
    
*  USEFUL EXPRESSION
      DRSQ=D*RSQ
*   
*  DEFAULT TO NO DISTORTION IF ILLEGAL KDIRN
      RFACTR=1D0
    
*  UNDISTORTED TO DISTORTED 
      IF(KDIRN.EQ.'UD') RFACTR=1D0+DRSQ 
    
*  DISTORTED TO UNDISTORTED 
      IF(KDIRN.EQ.'DU') RFACTR=(2D0*DRSQ+1D0)/(3D0*DRSQ+1D0)
    
*  ADJUST X, Y  
      X=X*RFACTR
      Y=Y*RFACTR
    
      END   
      SUBROUTINE SCHMIT(A,B,CURV)   
    
*  CALCULATE CURVATURE CORRECTION FOR SCHMIDT PLATES
*  X = XI / CURV ,   Y = ETA / CURV 
*  XI = X * CURV ,   ETA = Y * CURV 
    
    
*  GIVEN:   
*     A       TANGENTIAL COORD (XI) OR IDEAL PLATE COORD (X)
*     B       TANGENTIAL COORD (ETA) OR IDEAL PLATE COORD (Y)   
    
*  RETURNED:
*     CURV    CURVATURE CORRECTION FACTOR   
    
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION A,B,CURV 
    
    
      R=SQRT(A*A + B*B) 
      CURV=1D0  
      IF(ABS(R).GT.1D-20) CURV=(SIN(R)/COS(R))/R
    
      END   
      SUBROUTINE RMS(LUR,LUS,SIG,NREF)  
    
*  CALCULATE AND PRINTS RMS OF RESIDUALS
*  DR, DD, EW, SKY  
    
    
*  GIVEN:   
*     LUR     LOGICAL UNIT FOR FULL REPORT  
*     LUS     LOGICAL UNIT FOR SYNOPSIS 
*     NREF    NUMBER OF REFERENCE STARS 
    
*  MODIFIED:
*     SIG     RMS SUMMATION VECTOR  
    
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      INTEGER LUR,LUS,NREF  
      DOUBLE PRECISION SIG(4)   
    
*  ARC SECONDS TO RADIANS   
      PARAMETER (AS2R=4.8481368110953599D-6)
    
    
      REFN=DBLE(NREF)   
      DO I=1,4  
         SIG(I)=SQRT(SIG(I)/REFN)/AS2R  
      END DO
      WRITE(LUR,1000) SIG   
 1000 FORMAT(/75X,'RMS :',F9.2,2F10.2,F9.2/)
      WRITE(LUS,1001) SIG(3),SIG(2),SIG(4)  
 1001 FORMAT(/17X,'RMS :',F9.2,2F10.2)  
    
      END   
      SUBROUTINE DFLTIN(STRING,NSTRT,LENFLD,LENSTR,DRESLT,JFLAG)
*   
*   
*  CONVERT FREE-FORMAT INPUT INTO DOUBLE PRECISION FLOATING POINT   
*   
*  GIVEN:   
*     STRING      CHARACTER STRING CONTAINING FIELD TO BE DECODED   
*     NSTRT       POINTER TO 1ST CHARACTER OF FIELD IN STRING   
*     LENFLD      FIELD LENGTH (CHARACTERS) 
*     LENSTR      STRING LENGTH (CHARACTERS)
*   
*  RETURNED:
*     NSTRT       INCREMENTED   
*     DRESLT      DOUBLE PRECISION RESULT   
*     JFLAG       -1 = -OK, 0 = +OK, 1 = NULL FIELD, 2 = ERROR  
*   
*   
*     NOTES:
*        1     THE BASIC FORMAT IS #^.^@#^ WHERE # MEANS + OR -,
*               ^ MEANS A DECIMAL SUBFIELD AND @ MEANS D OR E.  
*        2     BLANKS:  
*               LEADING BLANKS ARE IGNORED. 
*               EMBEDDED BLANKS ARE ALLOWED ONLY AFTER # AND D OR E, AND
*               AFTER . WHERE THE FIRST ^ IS ABSENT.
*               TRAILING BLANKS ARE IGNORED; THE FIRST SIGNIFIES
*               END OF DECODING AND SUBSEQUENT ONES ARE SKIPPED.
*        3     FIELD SEPARATORS:
*               ANY CHARACTER OTHER THAN +,-,0-9,.,D,E OR BLANK MAY BE  
*               USED TO END A FIELD. COMMA IS RECOGNISED BY DFLTIN  
*               AS A SPECIAL CASE; IT IS SKIPPED, LEAVING THE   
*               POINTER ON THE NEXT CHARACTER. SEE 11, BELOW.   
*        4     BOTH SIGNS ARE OPTIONAL. THE DEFAULT IS +.   
*        5     THE MANTISSA DEFAULTS TO 1.  
*        6     THE EXPONENT DEFAULTS TO D0. 
*        7     THE DECIMAL SUBFIELDS MAY BE OF ANY LENGTH.  
*        8     THE DECIMAL POINT IS OPTIONAL FOR INTEGERS.  
*        9     A NULL FIELD IS ONE THAT DOES NOT BEGIN WITH 
*               +,-,0-9,.,D OR E, OR CONSISTS ENTIRELY OF BLANKS.   
*               IF THE FIELD IS NULL, JFLAG IS SET TO 1 AND 
*               RESULT IS LEFT UNTOUCHED.   
*        10    NSTRT = 1 FOR THE FIRST CHARACTER IN THE STRING. 
*        11    ON RETURN FROM DFLTIN, NSTRT POINTS TO THE CHARACTER 
*               FOLLOWING THE FINAL CHARACTER OF THE FIELD. IF A
*               SEPARATOR IS BEING USED - OTHER THAN COMMA - NSTRT  
*               MUST BE INCREMENTED BEFORE THE NEXT CALL TO DFLTIN. 
*        12    ERRORS (JFLAG=2) OCCUR WHEN: 
*              A) A +, -, D OR E IS LEFT UNSATISFIED.   
*              B) THE DECIMAL POINT IS PRESENT WITHOUT AT LEAST 
*                  ONE DECIMAL SUBFIELD.
*              C) AN EXPONENT MORE THAN 100 HAS BEEN PRESENTED. 
*               VARIOUS ARITHMETIC FAULTS ARE POSSIBLE WHICH WILL   
*               BE REPORTED BY THE ENVIRONMENT  
*               RATHER THAN BY DFLTIN.  
*        13    WHEN AN ERROR HAS BEEN DETECTED, NSTRT IS LEFT   
*               POINTING TO THE CHARACTER FOLLOWING THE LAST
*               ONE USED BEFORE THE ERROR CAME TO LIGHT. THIS   
*               MAY BE AFTER THE POINT AT WHICH A MORE SOPHISTICATED
*               PROGRAM COULD HAVE DETECTED THE ERROR. FOR EXAMPLE, 
*               DFLTIN DOES NOT DETECT THAT '1D999' IS UNACCEPTABLE 
*               UNTIL THE WHOLE FIELD HAS BEEN READ.
*        14    CERTAIN HIGHLY UNLIKELY COMBINATIONS OF MANTISSA &   
*               EXPONENT CAN CAUSE ARITHMETIC FAULTS DURING THE 
*               DECODE, DESPITE THE FACT THAT THEY TOGETHER 
*               COULD BE CONSTRUED AS A VALID NUMBER.   
*        15    DECODING IS LEFT TO RIGHT, ONE PASS. 
*        16    END OF FIELD MAY OCCUR IN ONE OF THREE WAYS: 
*              A) AS DICTATED BY THE FIELD LENGTH, LENFLD.  
*              B) AS DICTATED BY THE STRING LENGTH, LENSTR. 
*              C) DETECTED DURING THE DECODE.   
*              (C OVERRIDES B, AND B OVERRIDES A.)  
    
      IMPLICIT DOUBLE PRECISION (D) 
    
      CHARACTER*(*) STRING  
      INTEGER NSTRT,LENFLD,LENSTR   
      DOUBLE PRECISION DRESLT   
      INTEGER JFLAG 
    
      DATA DTEN,DTEN10/1D1,1D10/
    
    
    
    
    
    
    
*  SET UP POINTERS: 
*   CURRENT CHARACTER   
      JPTR=NSTRT
*   END OF FIELD (OR STRING)
      JNXFLD=MIN0(JPTR+LENSTR-1,LENFLD) 
    
*  SET DEFAULTS: MANTISSA & SIGN, EXPONENT & SIGN, DECIMAL PLACE COUNT  
      DMANT=0D0 
      MSIGN=1   
      NEXP=0
      ISIGNX=1  
      NDP=0 
    
*  LOOK FOR SIGN
 100  CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO ( 400,  100,  800,  500,  300,  200, 9110, 9100, 9110),NVEC 
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  NEGATIVE 
 200  CONTINUE  
      MSIGN=-1  
    
*  LOOK FOR FIRST LEADING DECIMAL   
 300  CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO ( 400, 300,  800,  500, 9200, 9200, 9200, 9200, 9210),NVEC  
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  ACCEPT LEADING DECIMALS  
 400  CONTINUE  
      DMANT=DMANT*DTEN+DIGIT
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO ( 400, 1310,  900,  600, 1300, 1300, 1300, 1300, 1310),NVEC 
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  LOOK FOR DECIMAL WHEN NONE PRECEDED THE POINT
 500  CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO ( 700, 500, 9200, 9200, 9200, 9200, 9200, 9200, 9210),NVEC  
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  LOOK FOR TRAILING DECIMALS   
 600  CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO ( 700, 1310,  900, 1300, 1300, 1300, 1300, 1300, 1310),NVEC 
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  ACCEPT TRAILING DECIMALS 
 700  CONTINUE  
      NDP=NDP+1 
      DMANT=DMANT*DTEN+DIGIT
      GO TO 600 
    
*  EXPONENT SYMBOL FIRST IN FIELD: DEFAULT MANTISSA TO 1
 800  CONTINUE  
      DMANT=1D0 
    
*  LOOK FOR SIGN OF EXPONENT
 900  CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO (1200, 900, 9200, 9200, 1100, 1000, 9200, 9200, 9210),NVEC  
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  EXPONENT NEGATIVE
 1000 CONTINUE  
      ISIGNX=-1 
    
*  LOOK FOR FIRST DIGIT OF EXPONENT 
 1100 CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO (1200, 1100, 9200, 9200, 9200, 9200, 9200, 9200, 9210),NVEC 
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  USE EXPONENT DIGIT   
 1200 CONTINUE  
      NEXP=NEXP*10+NDIGIT   
      IF(NEXP.GT.100) GO TO 9200
    
*  LOOK FOR SUBSEQUENT DIGITS OF EXPONENT   
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO (1200, 1310, 1310, 1310, 1310, 1310, 1310, 1310, 1300),NVEC 
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
*  COMBINE EXPONENT AND DECIMAL PLACE COUNT 
 1300 CONTINUE  
      JPTR=JPTR-1   
 1310 CONTINUE  
      NEXP=NEXP*ISIGNX-NDP  
    
*  SKIP IF NET EXPONENT NEGATIVE
      IF(NEXP.LT.0) GO TO 1500  
    
*  POSITIVE EXPONENT: SCALE UP  
 1400 CONTINUE  
      IF(NEXP.LT.10) GO TO 1410 
      DMANT=DMANT*DTEN10
      NEXP=NEXP-10  
      GO TO 1400
 1410 CONTINUE  
      IF(NEXP.LT.1) GO TO 1600  
      DMANT=DMANT*DTEN  
      NEXP=NEXP-1   
      GO TO 1410
    
*  NEGATIVE EXPONENT: SCALE DOWN
 1500 CONTINUE  
      IF(NEXP.GT.-10) GO TO 1510
      DMANT=DMANT/DTEN10
      NEXP=NEXP+10  
      GO TO 1500
 1510 CONTINUE  
      IF(NEXP.GT.-1) GO TO 1600 
      DMANT=DMANT/DTEN  
      NEXP=NEXP+1   
      GO TO 1510
    
*  GET RESULT & STATUS  
 1600 CONTINUE  
      J=0   
      IF(MSIGN.EQ.1) GO TO 1610 
      J=-1  
      DMANT=-DMANT  
 1610 CONTINUE  
      DRESLT=DMANT  
    
*  SKIP TO END OF FIELD 
 1620 CONTINUE  
      CALL IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)  
      GO TO (1720, 1620, 1720, 1720, 1720, 1720, 9900, 1720, 9900),NVEC 
*             0-9   SP    D/E    .     +     -     ,   ELSE   END   
    
 1720 CONTINUE  
      JPTR=JPTR-1   
      GO TO 9900
    
    
*  EXITS
    
*  NULL FIELD   
 9100 CONTINUE  
      JPTR=JPTR-1   
 9110 CONTINUE  
      J=1   
      GO TO 9900
    
*  ERRORS   
 9200 CONTINUE  
      JPTR=JPTR-1   
 9210 CONTINUE  
      J=2   
    
*  RETURN   
 9900 CONTINUE  
      NSTRT=JPTR
      JFLAG=J   
    
      END   
      SUBROUTINE IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT)
    
*   
*  IDENTIFY NEXT CHARACTER IN STRING
*   
*  GIVEN:   
*     STRING      CHARACTER STRING  
*     JPTR        POINTER TO CHARACTER TO BE IDENTIFIED 
*     JNXFLD      POINTER TO END OF FIELD   
*   
*  RETURNED:
*     JPTR        INCREMENTED UNLESS END OF FILED   
*     NVEC        VECTOR FOR IDENTIFIED CHARACTER   
*     NDIGIT      INTEGER 0-9 IF CHARACTER WAS A NUMERAL
*     DIGIT       DOUBLE PRECISION EQUIVALENT OF NDIGIT 
*   
*  NVEC TAKES THE FOLLOWING VALUES: 
*     1     0-9 
*     2     SPACE   
*     3     D,E 
*     4     .   
*     5     +   
*     6     -   
*     7     ,   
*     8     ELSE
*     9     END OF FIELD
    
    
      CHARACTER*(*) STRING  
      INTEGER JPTR,JNXFLD,NVEC,NDIGIT   
      DOUBLE PRECISION DIGIT
    
*  CHARACTER/VECTOR TABLES  
      CHARACTER KCTAB(17)   
      INTEGER KVTAB(17) 
      DATA KCTAB/'0','1','2','3','4','5','6','7','8','9',   
     *           ' ','D','E','.','+','-',','/   
      DATA KVTAB/10*1,2,2*3,4,5,6,7/
    
    
*  HANDLE END OF FIELD  
      IF(JPTR.GT.JNXFLD) THEN   
         NVEC=9 
      ELSE  
    
*     NOT END OF FIELD  
    
*     IDENTIFY CHARACTER
         DO NCHAR=1,17  
            IF(KCTAB(NCHAR).EQ.STRING(JPTR:JPTR)) GO TO 2200
         END DO 
    
*     UNRECOGNISED  
         NVEC=8 
         GO TO 2300 
    
*     RECOGNISED
 2200    CONTINUE   
         NVEC=KVTAB(NCHAR)  
    
*     ALLOW FOR NUMERALS
         NDIGIT=NCHAR-1 
         DIGIT=DBLE(FLOAT(NDIGIT))  
    
*     INCREMENT POINTER 
 2300    CONTINUE   
         JPTR=JPTR+1
      END IF
    
*  RETURN   
    
      END   
      SUBROUTINE DMATRX(N,A,Y,D,JF,W,IW)
*   
*  MATRIX INVERSION & SOLUTION TO SIMULTANEOUS EQUATIONS
*   
*  DOUBLE PRECISION 
*   
*  FOR THE SET OF N SIMULTANEOUS EQUATIONS IN N UNKNOWNS:   
*     A.Y = X   
*  WHERE:   
*     A IS A NON-SINGULAR N X N MATRIX  
*     Y IS THE VECTOR OF N UNKNOWNS 
*     X IS THE KNOWN VECTOR 
*  DMATRX COMPUTES: 
*     THE INVERSE OF MATRIX A   
*     THE DETERMINANT OF MATRIX A   
*     THE VECTOR OF N UNKNOWNS  
*   
*  ARGUMENTS:   
*     SYMBOL  TYPE DIMENSION           BEFORE              AFTER
*       N      I                   NO. OF UNKNOWNS       UNCHANGED  
*       A      DP    (N,N)             MATRIX             INVERSE   
*       Y      DP     (N)              VECTOR            SOLUTION   
*       D      DP                         -             DETERMINANT 
*     * JF     I                          -           SINGULARITY FLAG  
*       W      DP    (N,N)                -              WORKSPACE  
*       IW     I      (N)                 -              WORKSPACE  
*   
*  *  JF IS THE SINGULARITY FLAG.  IF THE MATRIX IS NON-SINGULAR,   
*    JF=0 IS RETURNED.  IF THE MATRIX IS SINGULAR, JF=-1 & D=0.0 ARE
*    RETURNED.  IN THE LATTER CASE, THE CONTENTS OF ARRAY A ON RETURN   
*    ARE UNDEFINED. 
*   
*  PEDIGREE:
*     BY PAUL MURDIN OUT OF UNIVERSITY OF ROCHESTER, NEW YORK.  
*   
*  ALGORITHM:   
*     NO DESCRIPTION AVAILABLE.  SEEMS TO BE BASED ON A STANDARD
*     METHOD BUT WITH RIGOUR IN THE SELECTION OF PIVOT POINTS   
*     COMPROMISED TO IMPROVE SPEED AND SIZE AT THE EXPENSE OF ACCURACY. 
*   
*  SPEED:   
*     VERY FAST 
*   
*  ACCURACY:
*     FAIRLY ACCURATE - ERRORS 1 TO 4 TIMES THOSE OF ROUTINES OPTIMISED 
*     FOR ACCURACY. 
*   
    
    
    
      DOUBLE PRECISION A(N,N),Y(N),D,W(N,N) 
      INTEGER IW(N) 
    
      DOUBLE PRECISION AMX,T,AKK,YK,AIK,SFA 
      PARAMETER (SFA=1D-20) 
    
    
      JF=0  
      D=1D0 
      DO K=1,N  
         AMX=DABS(A(K,K))   
         IMX=K  
         IF (K.NE.N) THEN   
            DO I=K+1,N  
               T=DABS(A(I,K))   
               IF (T.GT.AMX) THEN   
                  AMX=T 
                  IMX=I 
               END IF   
            END DO  
         END IF 
         IF (AMX.LT.SFA) THEN   
            JF=-1   
         ELSE   
            IF (IMX.NE.K) THEN  
               DO J=1,N 
                  T=A(K,J)  
                  A(K,J)=A(IMX,J)   
                  A(IMX,J)=T
               END DO   
               T=Y(K)   
               Y(K)=Y(IMX)  
               Y(IMX)=T 
               D=-D 
            END IF  
            IW(K)=IMX   
            AKK=A(K,K)  
            D=D*AKK 
            IF (DABS(D).LT.SFA) THEN
               JF=-1
            ELSE
               AKK=1D0/AKK  
               A(K,K)=AKK   
               DO J=1,N 
                  IF (J.NE.K) A(K,J)=A(K,J)*AKK 
               END DO   
               YK=Y(K)*AKK  
               Y(K)=YK  
               DO I=1,N 
                  AIK=A(I,K)
                  IF (I.NE.K) THEN  
                     DO J=1,N   
                        IF (J.NE.K) A(I,J)=A(I,J)-AIK*A(K,J)
                     END DO 
                     Y(I)=Y(I)-AIK*YK   
                  END IF
               END DO   
               DO I=1,N 
                  IF (I.NE.K) A(I,K)=-A(I,K)*AKK
               END DO   
            END IF  
         END IF 
      END DO
      IF (JF.NE.0) THEN 
         D=0D0  
      ELSE  
         DO K=1,N   
            NP1MK=N+1-K 
            KI=IW(NP1MK)
            IF (NP1MK.NE.KI) THEN   
               DO I=1,N 
                  T=A(I,NP1MK)  
                  A(I,NP1MK)=A(I,KI)
                  A(I,KI)=T 
               END DO   
            END IF  
         END DO 
      END IF
      END   
      SUBROUTINE DR2AF(NDP,ANGLE,SIGN,IDMSF)
*   
*   
*  CONVERT ANGLE FOR OUTPUT 
*  ------------------------ 
*   
*   
*  GIVEN:   
*     NDP           NUMBER OF DECIMAL PLACES OF ARC SECONDS 
*     ANGLE         ANGLE IN RADIANS (DOUBLE PRECISION) 
*   
*  RETURNED:
*     SIGN          '+' OR '-'  
*     IDMSF         DEGREES, ARC MINUTES, ARC SECONDS, FRACTION 
*   
*   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
*   
      INTEGER NDP   
      DOUBLE PRECISION ANGLE
      CHARACTER*1 SIGN  
      INTEGER IDMSF(4)  
    
*  RADIANS TO ARC SECONDS   
      DOUBLE PRECISION R2AS 
      PARAMETER (R2AS=206264.80624709636D0) 
    
*  HANDLE SIGN  
      SIGN='+'  
      IF(ANGLE.LT.0D0) SIGN='-' 
    
*  FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE 
      RAS=REAL(10**MAX(NDP,0))  
      RAM=RAS*60D0  
      RD=RAM*60D0   
    
*  ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED   
      A=NINT(RAS*R2AS*ABS(ANGLE))   
    
*  SEPARATE INTO FIELDS 
      AD=INT(A/RD)  
      A=A-AD*RD 
      AM=INT(A/RAM) 
      A=A-AM*RAM
      AS=INT(A/RAS) 
      AF=A-AS*RAS   
    
*  RETURN RESULTS   
      IDMSF(1)=INT(AD)  
      IDMSF(2)=INT(AM)  
      IDMSF(3)=INT(AS)  
      IDMSF(4)=INT(AF)  
    
      END   
      SUBROUTINE DR2TF(NDP,ANGLE,SIGN,IHMSF)
*   
*   
*  CONVERT ANGLE FOR OUTPUT, EXPRESSED AS TIME  
*  -------------------------------------------  
*   
*   
*  GIVEN:   
*     NDP           NUMBER OF DECIMAL PLACES OF SECONDS 
*     ANGLE         ANGLE IN RADIANS (DOUBLE PRECISION) 
*   
*  RETURNED:
*     SIGN          '+' OR '-'  
*     IHMSF         HOURS, MINUTES, SECONDS, FRACTION   
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      INTEGER NDP   
      DOUBLE PRECISION ANGLE
      CHARACTER*1 SIGN  
      INTEGER IHMSF(4)  
    
*  RADIANS TO SECONDS   
      DOUBLE PRECISION R2S  
      PARAMETER (R2S=13750.987083139757D0)  
    
*  HANDLE SIGN  
      SIGN='+'  
      IF(ANGLE.LT.0D0) SIGN='-' 
    
*  FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE 
      RS=REAL(10**MAX(NDP,0))   
      RM=RS*60D0
      RH=RM*60D0
    
*  ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED   
      A=NINT(RS*R2S*ABS(ANGLE)) 
    
*  SEPARATE INTO FIELDS 
      AH=INT(A/RH)  
      A=A-AH*RH 
      AM=INT(A/RM)  
      A=A-AM*RM 
      AS=INT(A/RS)  
      AF=A-AS*RS
    
*  RETURN RESULTS   
      IHMSF(1)=INT(AH)  
      IHMSF(2)=INT(AM)  
      IHMSF(3)=INT(AS)  
      IHMSF(4)=INT(AF)  
    
      END   
      DOUBLE PRECISION FUNCTION DRANGE(ANGLE)   
*   
*   
*  NORMALISE ANGLE INTO RANGE +/- PI
*  ---------------------------------
*   
*   
*  ANGLE IS THE ANGLE IN RADIANS
*  BOTH ANGLE AND THE RESULT ARE DOUBLE PRECISION   
*   
    
      DOUBLE PRECISION ANGLE
    
      DOUBLE PRECISION DPI,D2PI 
      PARAMETER (DPI=3.141592653589793238462643D0)  
      PARAMETER (D2PI=6.283185307179586476925287D0) 
    
    
      DRANGE=DMOD(ANGLE,D2PI)   
      IF(DABS(DRANGE).GE.DPI) DRANGE=DRANGE-DSIGN(D2PI,ANGLE)   
    
      END   
      DOUBLE PRECISION FUNCTION DRANRM(ANGLE)   
*   
*   
*  NORMALISE ANGLE INTO RANGE 0-2PI 
*  -------------------------------- 
*   
*   
*  ANGLE IS THE ANGLE IN RADIANS
*  BOTH ANGLE AND THE RESULT ARE DOUBLE PRECISION   
    
      DOUBLE PRECISION ANGLE
    
      DOUBLE PRECISION D2PI 
      PARAMETER (D2PI=6.283185307179586476925287D0) 
    
    
      DRANRM=DMOD(ANGLE,D2PI)   
      IF(DRANRM.LT.0D0) DRANRM=DRANRM+D2PI  
    
      END   
      SUBROUTINE DS2TP(RA,DEC,RAZ,DECZ,XI,ETA)  
*   
*   
*  PROJECTION OF SPHERICAL COORDINATES ONTO TANGENT PLANE   
*  ('GNOMONIC' PROJECTION - 'STANDARD COORDINATES') 
*   
*   
*  GIVEN:   
*     RA,DEC      SPHERICAL COORDINATES OF POINT TO BE PROJECTED
*     RAZ,DECZ    SPHERICAL COORDINATES OF TANGENT POINT
*   
*  RETURNED:
*     XI,ETA      RECTANGULAR COORDINATES ON TANGENT PLANE  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION RA,DEC,RAZ,DECZ,XI,ETA   
    
    
      SDECZ=SIN(DECZ)   
      SDEC=SIN(DEC) 
      CDECZ=COS(DECZ)   
      CDEC=COS(DEC) 
      RADIF=RA-RAZ  
      SRADIF=SIN(RADIF) 
      CRADIF=COS(RADIF) 
    
      DENOM=SDEC*SDECZ+CDEC*CDECZ*CRADIF
    
      XI=CDEC*SRADIF/DENOM  
      ETA=(SDEC*CDECZ-CDEC*SDECZ*CRADIF)/DENOM  
    
      END   
      SUBROUTINE DTP2S(XI,ETA,RAZ,DECZ,RA,DEC)  
*   
*   
*  TRANSFORM TANGENT PLANE COORDINATES INTO SPHERICAL   
*   
*   
*  GIVEN:   
*     XI,ETA      TANGENT PLANE RECTANGULAR COORDINATES 
*     RAZ,DECZ    SPHERICAL COORDINATES OF TANGENT POINT
*   
*  RETURNED:
*     RA,DEC      SPHERICAL COORDINATES 
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION XI,ETA,RAZ,DECZ,RA,DEC   
    
    
      SDECZ=SIN(DECZ)   
      CDECZ=COS(DECZ)   
    
      DENOM=CDECZ-ETA*SDECZ 
      RADIF=ATAN2(XI,DENOM) 
    
      RA=DRANRM(RADIF+RAZ)  
      DEC=ATAN2(SDECZ+ETA*CDECZ,SQRT(XI*XI+DENOM*DENOM))
    
      END   
      SUBROUTINE MMP(RM0, DM0, PR, PD, EPOCH0, EPOCH1, RM1, DM1)
*   
*  GIVEN:   
*     RM0, DM0           MEAN RA, DEC AT 1ST EPOCH  
*     PR, PD             PROPER MOTIONS IN RA, DEC AT SPECIFIED EPOCH   
*                          (RADIANS PER YEAR)   
*     EPOCH0             1ST EPOCH  
*     EPOCH1             2ND EPOCH  
*   
*  RETURNED:
*     RM1, DM1           MEAN RA, DEC AT 2ND EPOCH  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION RM0, DM0, PR, PD, EPOCH0, EPOCH1, RM1, DM1   
    
      DOUBLE PRECISION RMATP(3,3), V(3) 
    
    
*  PROPER MOTIONS   
    
      TAU= EPOCH1 - EPOCH0  
      R= RM0 + PR*TAU   
      D= DM0 + PD*TAU   
    
    
*  PRECESSION MATRIX
    
      CALL PREC(EPOCH0, EPOCH1, RMATP)  
    
    
*  MEAN TO MEAN 
    
      CALL DCS2C(R, D, V)   
      CALL RMATMU(V, RMATP) 
      CALL DCC2S(R, DM1, V) 
      RM1= DRANRM(R)
    
      END   
      SUBROUTINE PREC(EPOCH0, EPOCH1, S)
*   
*  GIVEN:   
*     EPOCH0             STARTING EPOCH (YEARS AD)  
*     EPOCH1             ENDING EPOCH (YEARS AD)
*   
*  RETURNED:
*     S                  PRECESSION MATRIX  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION EPOCH0, EPOCH1, S(3,3)   
    
    
      TSER(C0,C1,C2)=DRANGE(C0+(C1+C2*T)*T) 
    
    
*  TIME FROM 1900.0 TO ORIGINAL EPOCH (TROPICAL CENTURIES)  
    
      T0=(EPOCH0 - 1900D0)/100D0
    
    
*  TIME SINCE ORIGINAL EPOCH (TROPICAL CENTURIES)   
    
      T= (EPOCH1 - EPOCH0)/100D0
      TSQ=T*T   
    
    
*  PRECESSION PARAMETERS
    
      Z0=TSER(0., 6.768D-6*T0 + 1.11713192D-2, 1.464D-6)
      Z=Z0 + 3.835D-6*TSQ   
      TH=TSER(0D0, -4.135D-6*T0 + 9.7189726D-3, -2.065D-6)  
    
      SZ0=SIN(Z0)   
      CZ0=COS(Z0)   
    
      SZ=SIN(Z) 
      CZ=COS(Z) 
    
      STH=SIN(TH)   
      CTH=COS(TH)   
    
      CTHSZ=CTH*SZ  
      CTHCZ=CTH*CZ  
    
    
*  ROTATION MATRIX  
    
      S(1,1)= CZ0*CTHCZ - SZ0*SZ
      S(1,2)=-SZ0*CTHCZ - CZ0*SZ
      S(1,3)=-STH*CZ
      S(2,1)= CZ0*CTHSZ + SZ0*CZ
      S(2,2)=-SZ0*CTHSZ + CZ0*CZ
      S(2,3)=-STH*SZ
      S(3,1)= CZ0*STH   
      S(3,2)=-SZ0*STH   
      S(3,3)= CTH   
    
      END   
      SUBROUTINE PMC(PR, PD, EPOCH, DATE, DR, DD)   
*   
*  GIVEN:   
*     PR, PD             PROPER MOTIONS IN RA, DEC AT SPECIFIED EPOCH   
*                          (RADIANS PER YEAR)   
*     EPOCH              YEAR AD OF MEAN EQUATOR & EQUINOX  
*     DATE               DAYS SINCE JD 2400000.0 FOR APPARENT PLACE 
*   
*  RETURNED:
*     DR, DD             TOTAL PROPER MOTION CORRECTIONS (RADIANS)  
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION PR, PD, EPOCH, DATE, DR, DD  
    
    
*  TIME SINCE ORIGINAL EPOCH (TROPICAL YEARS)   
    
      TAU= (DATE - 15020.313D0)/365.2421957D0 - (EPOCH - 1900D0)
    
    
*  PROPER MOTION CORRECTIONS (RA,DEC VECTOR)
    
      DR= PR*TAU
      DD= PD*TAU
    
      END   
      SUBROUTINE DCS2C(A,B,V)   
*   
*     A IS RA/HA/LONG/AZ
*     B IS DEC/LAT/ZD   
*     V IS X,Y,Z
*   
    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION A, B, V(3)   
    
    
      AR= DRANGE(A) 
      BR= DRANGE(B) 
    
      COSBR= COS(BR)
    
      V(1)= COS(AR)*COSBR   
      V(2)= SIN(AR)*COSBR   
      V(3)= SIN(BR) 
    
      END   
      SUBROUTINE DCC2S(A,B,V)   
*   
*     A IS RA/HA/LONG/AZ
*     B IS DEC/LAT/ZD   
*     V IS X,Y,Z
*   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
    
      DOUBLE PRECISION A, B, V(3)   
    
    
      X=V(1)
      Y=V(2)
    
      A=ATAN2(Y,X)  
      B=ATAN2(V(3), SQRT(X*X + Y*Y))
    
      END   
      SUBROUTINE RMATMU(A,B)
*   
*  PERFORMS THE 3-D FORWARD UNITARY TRANSFORMATION: 
*     VECTOR A = MATRIX B * VECTOR A
*   
    
      DOUBLE PRECISION A(3), B(3,3) 
    
      DOUBLE PRECISION C(3) 
    
    
      DO J=1,3  
         C(J)=0D0   
         DO I=1,3   
            C(J)=C(J)+A(I)*B(J,I)   
         END DO 
      END DO
    
      DO J=1,3  
         A(J)=C(J)  
      END DO
    
      END   
      FUNCTION RJDATE ( NYEAR, NDIY )   
*   
*   
*  JULIAN DATE FROM YEAR & DAY-IN-YEAR  
*   
*   
*  ARGUMENTS:   
*   
*        NYEAR GREGORIAN YEAR AD (INTEGER; MUST BE 1901 THROUGH 2100)   
*        NDIY  DAY IN YEAR NUMBER (INTEGER; JAN 1ST = DAY 1)
*   
*        RESULT IS THE JULIAN DATE MINUS 2400000.0 FOR 0 HRS UT 
*   
    
      RJDATE=REAL(NDIY+(NYEAR-1)/4)+REAL(NYEAR-1901)*365.0+14909.5  
    
      END   


