C Current Version: Jan 17, 2000
C "cuepoch" test for 2000 changed!
C
C whydra.f
C
C Phil's Pretty Good Assignment Routine
C written in Archaic FORTRAN
C based in part upon Tod Lauer's assignment algorithm
C with Taft's changes
C but provided with Returns to the Well.
C Begun this day 28 April 1991
C Zeroth order version finished May 22, 1991.
C Begans mods to 12/19/93 version for WIYN; started May 6, 1994!
C First compiled test version May 18, 1994.
C Major mod: added "refit" routine May 25, 1994
C
C Major mods 17 March 1995:
C a) incorporate Lockheed astrometry and PLM modeling
C b) correct problem with rotator convention in "getangle"
C c) added Sam's new warp-plate values
C
C added the "rest of the astrometry" on 24 March 1995 at NAU.
C modified 8/28/95 to (a) make plate scale 9.400 to conform with Sam
C set unit=9 to scratch
C added red fibers 276 and 279
C
C BUG FIX 9/18/95: now will only update to new center and rotation angle
C if the number of required fops and skys are met.  Also added a warning
C statement to the assignment output if this is not the case.
C 
C I also got rid of that annoying 'HUH?' and replaced it by 'EDGE'
C
C 9/26/95 Commented out all writes to unit 9.  Program now does not
C require so much disk space.
C 
C I have also removed the ancient reference to "acquisition stars".
C
C changes 2 October
C changed "park circle" radius to 76200 from 76300.  Although the park
C circle is at 76300, no buttons can be assigned further than 76200
C encoder units.
C this was the true source of the "HUH?" statment above.
C 
C Also 2 October:
C redefined "setfibphys" to "setfibphysphil" to set to physical location
C  plus or minus 1.
C
C BUG FIX: 25 oct---added "Nfops" to itotal.
C 
C
C NEW BUG FIX 8 August 1997: problem in "gerard" index discovered by R. Wolff
C also added "nwant" before first call to "wings". 
C NEW BUG FIX 10 August 1997, again thanks to RW: call park(but,0)--->(ibut,0)
C down in "finalverify".  Also improved the efficiency of "legal".
C
C Modified "legal" again as per RW's code so it's a bit more efficient.
C ---Aug 12.  
C modified latitudes so they all agree----9/15/97
C several more minor bug fixes re RW---Oct 23, 1997 
C
C REPLACED METHOD FOR DETERMINING GOOD/BAD FIBERS Nov 29, 1999
C Now reads from concenticities file.....
C Unless it was tweaked! Stupid bug fixed 27 Feb 2000
	implicit real*8 (a-h,o-z)
	real*4 aweight
C data files
	character*1 ians
	character*1 itype
	character*4 stat
	character*80 filein,hydraout,logfile
	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
C it is also the maximum number of acquisition stars.
	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, fops, and acq.
C (The "acq" is actually an unused catagory, but may become parasitic
C fibers some day)
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),
     .      raacqs(maxguides),decacqs(maxguides),epochacqs(maxguides)
	dimension idobjs(maxstars),idskys(maxstars),idfops(maxguides),
     .            idacqs(maxguides),
     .            ifibs(maxfibers),ifopsfibs(maxguides)
        character*20 nameobjs(maxstars),nameskys(maxstars),
     .               nameacqs(maxguides),namefops(maxguides),namecent
	character*64 fieldname
	integer IDMSt(4),IDMSa(4)
	integer*4 objstatus(maxstars),skystatus(maxstars),
     .            acqstatus(maxguides),fopstatus(maxguides),
     .            xhop(maxfibers),yhop(maxfibers),hopstatus(maxfibers),
     .            fibstatus(maxfibers),inum(maxfibers),
     .            ixobjs(maxstars),iyobjs(maxstars),
     .            ixfops(maxguides),iyfops(maxguides),
     .	          ixskys(maxstars),iyskys(maxstars),
     .            ixacqs(maxguides),iyacqs(maxguides),
     .            fopsfibstatus(maxfibers),
     .		  itrackfops(maxguides),itrackacqs(maxguides),
     .		  itrackskys(maxstars),itrackobjs(maxstars)
C
C
C
C initiate the fiber file:
	call libinit
C
C
C
C This program takes an input file ("filein") of coordinates in either
C Old Nessie Format ---or----New Fangled Hydra Format (NFHF)
C     ---unit=1------
C
C
C
C ....and outputs the following 3 files:
C
C
C
C 
C
C "permanent" output file:
C    ---unit=2----
C              NFHF file containing the assignments and reasons for non
C                   assignements (all objects retained).  The ordering of
C                   the objects have been changed so that the assigned 
C  		    objects have been rotated to the bottom of the list.
C   ----unit=8----
C             "assignmentlog" which keeps a running account of all your
C              runs.  This file gets appended to every time you run
C              the program.
C   ----unit=9----
C               LOG file containing the very nitty-gritty issues of the
C               the assignments.  Mainly useful for debugging purposes.
C               These write statements have now been commented out.
C
CC Open input data file
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
	hydraout=out(1:i-1)//'.hydra'
	logfile=out(1:i-1)//'.log'
	open(unit=2,file=hydraout,status='new',err=5)
	open(unit=8,file='assignmentlog',access='append')
 	open(unit=9,status='scratch')
	write(*,*)hydraout
	write(8,*)hydraout
	go to 15
5	write(*,6)
	write(8,6)
6	format('Im sorry, but I cant open the output 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 filetype(ifile)
C ifile=0 Nessie file
C ifile=1 Hydra file
	Call grabvitals(fieldname,deepoch,cuepoch,aLST,exptime,
     .wave,camwave,icable,aweight,ifile)
C
C
C
C
C 
C
C Grab the RAs and DECs of all objects, classifying them at the
C same time
20	continue
	call grabradec(racent,deccent,epochcent,idcent,namecent,
     .  rafops,decfops,epochfops,idfops,namefops,Nfops,itrackfops,
     .   maxguides,
     .  raacqs,decacqs,epochacqs,idacqs,nameacqs,Nacqs,itrackacqs,
     .   maxguides,
     .  raskys,decskys,epochskys,idskys,nameskys,Nskys,itrackskys,
     .   maxstars,
     .  raobjs,decobjs,epochobjs,idobjs,nameobjs,Nobjs,itrackobjs,
     .   maxstars,
     .                 deepoch,ifile)
C bug here corrected!
	itotal=1+Nfops+Nacqs+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
Cwrite(*,21) Nacqs
Cwrite(8,21) Nacqs
21	format(1x,'There are',i5,' acquisition stars found')
C
C
	write(*,23) Nfops
	write(8,23) Nfops
23	format(1x,'There are',i5,' fops stars found')
	write(*,24)
	write(8,24)
24	format(1x,'How many fops need to be assigned to fops stars?: ',$)
	read(*,*) nfopsrequired
	write(8,*) nfopsrequired
	if(nfopsrequired.gt.Nfops) nfopsrequired=Nfops
	if(nfopsrequired.lt.0) nfopsrequired=0
C
C
	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')
	write(*,28)
	write(8,28)
28	format(1x,'How many skies are required?: ',$)
	read(*,*) nskysrequired
	write(8,*)nskysrequired
	if(nskysrequired.gt.Nskys) nskysrequired=Nskys
	if(nskysrequired.lt.0) nskysrequired=0
C determine the angle of N relative to fiber 0 when the telescope
C is on the meridian and at physical index 0.
C first precess to current epoch:
	call mmp(racent,deccent,0d0,0d0,epochcent,cuepoch,raz,dcz)
	call getangle(dcz,angle)
C ask the angle question
	iaweak=1
	write(*,651) angle
	write(8,651) angle
651	format(1x,'The default angle of N rel to fib 0 is',f8.3,' degs',
     &  /1x,'Would you like to tweak the angle? y or n: ',$)
	read(*,31) ians
	write(8,31) ians
	if((ians.eq.'y').or.(ians.eq.'Y')) then
644	write(*,652)
	write(8,652)
652	format(1x,'How many times? (3, 5, 7...) : ',$)
	read(*,*) iaweak
	write(8,*) iaweak
	if(iaweak.lt.1) iaweak=1
	if(iaweak.gt.15) then
	write(*,643)
	write(8,643)
643	format('Too many times: keep it to 15 or less or it will take',
     &  ' forever!')
	go to 644
	end if
	if( (2*(iaweak/2)+1) .ne. iaweak) then
	write(6,37)
	write(8,37)
	go to 644
	end if
667	write(*,665)
	write(8,665)
665	format('By how much? (Suggest 0.1 deg...): ',$)
	read(*,*) deltang
	write(8,*) deltang
	if((deltang.lt.0.01).or.(deltang.gt.20.)) then
	write(*,33)
	write(8,33)
	go to 667
	end if
	end if
C ask the center question
	itweak=1
29	write(*,30)
	write(8,30)
30	format(1x,'Would you like to tweak the center? y or n: ',$)
	read(*,31) ians
	write(8,31) ians
31	format(a1)
	if((ians.eq.'y').or.(ians.eq.'Y')) then
36	write(*,35)
	write(8,35)
35	format(1x,'How many grid points per side? (3, 5, 7...): ',$)
	read(*,*) itweak
	write(8,*) itweak
	if(itweak.lt.1) go to 29
	if( (2*(itweak/2)+1) .ne.itweak) then
	write(*,37)
	write(8,37)
37	format('Number needs to be odd')
	go to 36
	end if
	ngrid=itweak*itweak
38	write(*,32) ngrid
	write(8,32) ngrid
32	format(1x,'How many arcmin per each of the',i3,' grid points?: ',$)
	read(*,*) tweak
	write(8,*) tweak
		if((tweak.lt.0.01).or.(tweak.gt.60.)) then
			write(*,33)
			write(8,33)
33			format(1x,'Step size too big or small!')
			go to 38
		end if
	write(*,34) itweak,itweak,tweak
	write(8,34) itweak,itweak,tweak
	tweakrads=(3.14159/180.)*tweak/60.
34	format(1x,'Will check a',i3,'x',i3,' grid with step size',
     .f7.3,' arcmin')
	end if
C specify that the best centers are the original centers (true whether we loop
C or not)
	rabest=racent
	decbest=deccent
	angbest=angle
	ansbest=0.0d0
5000	continue
C We will now execute the following giant loop either once or ngrid times
C depending upon the value of itweak.
	icen=itweak/2+1
	iacen=iaweak/2+1
	do 5001 idec=1,itweak
	do 5001 ira=1,itweak
	do 5001 iangle=1,iaweak
	racentn=racent+(ira-icen)*tweakrads/dcos(deccent)	
	deccentn=deccent+(idec-icen)*tweakrads
	dangle=angle+(iangle-iacen)*deltang
Cwrite(9,5991) ira,racent,racentn,rabest,ansbest,
C    .               idec,deccent,deccentn,decbest,ansbest,
C    .               iangle,angle,dangle,angbest,ansbest
5991	format('in centering loop:'/10x,'ra:',i4,1x,3g12.6,f8.2,
     .                            /10x,'dc:',i4,1x,3g12.6,f8.2,
     .                           /10x,'an:',i4,1x,3f8.3,f8.2)
C
C
C
C
CC initialiate the fiber data-base and find out a few things:
C
	call initallfibers(icable,nfibs,ifibs,nfopsfibs,ifopsfibs)
C
C find out number of available fibers for this cable (nfibs)
C and fiber numbers (ifibs(i),i=1,nfibs)
C Plus the number of working fop-fibers (nfopsfibs) and their fiber numbers
C (ifopsfibs(i),i=1,nfopsfibs)

C
C

C
C
C Convert all objects to x,y "sky" coordinates, corrected for
C precession, refraction, distortion within the field-corrector,
C and the rotation angle.  Thus, the only difference beteen
C "sky" coordinates and "plate coorditinates" is flexure within
C  the gripper stages, which Lee implictly takes out.
C 
C
C
C
C
C
Cwrite(9,2681)
2681	format(1x,'Calling astrometry for OBJECTS')
	call wiynness(racentn,deccentn,epochcent,
     .               raobjs,decobjs,epochobjs,nobjs,ixobjs,iyobjs,
     .               objstatus,
     .               cuepoch,aLST,exptime,camwave,wave,
     .               dangle)
Cwrite(9,2682)
2682	format(1x,'Calling astrometry for SKYS')
	call wiynness(racentn,deccentn,epochcent,
     .               raskys,decskys,epochskys,nskys,ixskys,iyskys,
     .	             skystatus,
     .               cuepoch,aLST,exptime,camwave,wave,
     .               dangle)
Cwrite(9,2683)
2683	format(1x,'Calling astrometry for FOPS')
	call wiynness(racentn,deccentn,epochcent,
     .               rafops,decfops,epochfops,nfops,ixfops,iyfops,
     .               fopstatus,
     .               cuepoch,aLST,exptime,camwave,camwave,
     .               dangle)
Cwrite(9,2684)
2684	format(1x,'Calling astrometry for ACQS')
	call wiynness(racentn,deccentn,epochcent,
     .               raacqs,decacqs,epochacqs,nacqs,ixacqs,iyacqs,
     .               acqstatus,
     .               cuepoch,aLST,exptime,camwave,camwave,
     .               dangle)
	call dr2tf(3,racentn,sign,IDMSt)
	call dr2af(2,deccentn,sign,IDMSa)
	write(*,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle
	write(8,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle
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)
Cwrite(9,102) idobjs(i),nameobjs(i),IDMSt,sign,IDMSa,epochobjs(i),
C    .               ixobjs(i),iyobjs(i),objstatus(i)
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)
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)
Cwrite(9,102) idskys(i),nameskys(i),IDMSt,sign,IDMSa,epochskys(i),
C    .               ixskys(i),iyskys(i),skystatus(i)
1101	continue
1111	if(nfops.lt.1) go to 112
Cwrite(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)
Cwrite(9,102) idfops(i),namefops(i),IDMSt,sign,IDMSa,epochfops(i),
C    .                ixfops(i),iyfops(i),fopstatus(i)
105	continue
112	if(nacqs.lt.1) go to 113
Cwrite(9,107)
107	format(1x,'acqs:')
	do 108 i=1,nacqs
	call dr2tf(3,raacqs(i),sign,IDMSt)
	call dr2af(2,decacqs(i),sign,IDMSt)
Cwrite(9,102) idacqs(i),nameacqs(i),IDMSt,sign,IDMSa,epochacqs(i),
C    .                ixacqs(i),iyacqs(i),acqstatus(i)
108	continue
113	continue
C

C
C
C "I don't care if you did drive through a mountain this afternoon.  This
C is New Jersey."
C
C
C
C
C
C 
C
C
C Fill up Tod's hopper with the higest priority objects, AFTER checking
C to see if the position is "legal".
C
C
C
200	continue
C fill up Tod's hopper rather carefully.
C Sigh, I wish I didnt have to be so smart.
	if (nobjs.lt.1) go to 300
	ihop=0
	iobj=0
	nlast=0
202	continue
	if(ihop.eq.nfibs) go to 205
C stop when we are out of fibers
201	iobj=iobj+1
	if(iobj.gt.nobjs) go to 205
C stop when we are out of objects
	if(objstatus(iobj).lt.-1) go to 201
C don't put any illegal objects into the hopper
C
	ihop=ihop+1
	xhop(ihop)=ixobjs(iobj)
	yhop(ihop)=iyobjs(iobj)
	hopstatus(ihop)=-1
	inum(ihop)=iobj
C keep track of who we're putting in there!
	nlast=iobj
C keep track of the last one put in
	go to 202
205	Continue
C
C
C
C now do the preliminary assignment of the first ihop legal objects
C (ihop <= min0(nfibs,number-of-legal-objects)
	if(ihop.lt.1) go to 5001
	call tod(xhop,yhop,ihop,hopstatus,nfibs,ifibs,fibstatus)
C Whew!!!!!!!!!!!!!!!!!!!!
C 
C
C
C
C
C
C
C
300	continue
C
C Now set up the fops:
	if (nfops.lt.1) go to 400
	if (nfopsfibs.lt.1) go to 400
C load up the integer hopper having the x and y's of the potential fop stars
C initiate the status of the fops' fibers to -1
	do 302 i=1,nfopsfibs
	fopsfibstatus(i)=-1
302	continue
Cwrite(9,303) (ixfops(i),iyfops(i),fopstatus(i),i=1,nfops)
303	format(1x,'fop objects=',3i9)
Cwrite(9,304) (ifopsfibs(i),fopsfibstatus(i),i=1,nfopsfibs)
304	format(1x,'fop fibers=',2i9)
C Put the fops on for free if we can
C Change in case one wants ALL the fops assigned... 8/8/97---plm
	nwant=nfopsrequired
	if(nwant.gt.nfopsfibs-1) nwant=nfopsfibs-1
	call wings(ixfops,iyfops,nfops,fopstatus,
     .            nfopsfibs,ifopsfibs,fopsfibstatus,nwant)
C
C There is a potential problem here---could have ALL fops assigned.
C But only if they went on "for free".
C
C
C
C
C
C Now force the issue!
	call gerard(ixfops,iyfops,nfops,fopstatus,nfopsrequired,
     .nfopsfibs,ifopsfibs,fopsfibstatus,
     .nfibs,ifibs,fibstatus,hopstatus,
     .nfopsfibs,ifopsfibs,fopsfibstatus,aweight)
C
C
	call countem(ihop,hopstatus,nassigned,aweight,wts)
Cwrite(9,5004) nassigned,ihop
5004	format(1x,'After putting on fops we had',i5,' objs out of',i5,
     .     ' hopper objs')
C
C
400	continue
C
C
C
C Now do skies!
	if(nskys.lt.1) go to 500
	if(nfibs.lt.1) go to 500
Cwrite(9,405) 
405	format(1x,'Trying to sneak on some skys now:')
	call wings(ixskys,iyskys,nskys,skystatus,
     .             nfibs,ifibs,fibstatus,nskysrequired)
C
C
C
C
C Now force the issue:
Cwrite(9,406)
406	format(1x,'Trying to force the issues with skys now:')
	call gerard(ixskys,iyskys,nskys,skystatus,nskysrequired,
     .  nfibs,ifibs,fibstatus,
     .  nfibs,ifibs,fibstatus,hopstatus,
     .  nfopsfibs,ifopsfibs,fopsfibstatus,aweight)
C
C
	call countem(ihop,hopstatus,nassigned,aweight,wts)
Cwrite(9,5006) nassigned,ihop
5006	format(1x,'After putting on skys we had',i5,' objs out of',i5,
     .             ' hopper objs')

C Now that we are done fooling around with hopstatus, let's record the
C the actual objects as done:
500	continue
	if(nobjs.lt.1) go to 5001
	if(ihop.lt.1) go to  5001
	do 501 i=1,ihop
	objstatus(inum(i))=hopstatus(i)
501	continue
C Whew!

	call wings (ixobjs,iyobjs,nobjs,objstatus,
     .              nfibs,ifibs,fibstatus,nolimit)
C
C
	call countem(nobjs,objstatus,nassigned,aweight,wts)
Cwrite(9,5007) nassigned,nobjs
5007	format(1x,'After adding in objs we had',i5,' obs out of total',i5)
C
C Now try to put on any extra skys for free
	call wings(ixskys,iyskys,nskys,skystatus,
     .             nfibs,ifibs,fibstatus,nolimit)

C Now try to put on any extra fops for free
	call wings(ixfops,iyfops,nfops,fopstatus,
     .            nfopsfibs,ifopsfibs,fopsfibstatus,nolimit)



C finally, eliminate any possibility of a illegal fiber!
	call finalverify(nfibs,ifibs,fibstatus,objstatus,nobjs,
     .                   skystatus,nskys,
     .        nfopsfibs,ifopsfibs,fopsfibstatus,fopstatus,nfops)

	call countem(nfops,fopstatus,nf,aweight,wtsd)
	call countem(nskys,skystatus,ns,aweight,wtsd)
	call countem(nobjs,objstatus,nassigned,aweight,wts)
	write(*,5008) nassigned,nobjs,ns,nskys,nf,nfops,wts
	write(8,5008) nassigned,nobjs,ns,nskys,nf,nfops,wts
5008	format(1x,'Assigned: ' i4,'/',i4,' objs;'
     .i4,'/',i4,' skys;',i4,'/',i4,' fops;', '  Total wts=',f7.2)
C
C 
C Note if there aren't enough fops or skys:
	if(nf.lt.nfopsrequired) then
	write(*,5009)
	write(8,5009)
	end if
	if(ns.lt.nskysrequired) then
	write(*,5010)
	write(8,5010)
	end if
5009	format('Not enough FOPS assigned!')
5010	format('Not enough SKYS assigned!')

C
C
C
C
C

	if ((ira.eq.icen).and.(idec.eq.icen).and.(iangle.eq.iacen)) 
     &       anscent=wts
	if ((wts.gt.ansbest).and.(nf.ge.nfopsrequired).and.
     &     (ns.ge.nskysrequired)) then
		rabest=racentn
		decbest=deccentn
		angbest=dangle
		ansbest=wts
	end if
C Only update centers if (a) the weights are better, (b) the number of
C fops equals or exceeds the required number of fops, and (c) the number of skys
C equals or exceeds the number of skys required.
Cwrite(9,5992) ira,idec,iangle,wts,ansbest,rabest,decbest,angbest
5992	format(1x,'End loop: ira,idec,iangle,wts,ansbest,rabest,',
     .'decbest,angbest=',
     ./10x,3i5,2f8.2,2g12.6,f8.2)
5001	continue
	if((itweak.gt.1).or.(iaweak.gt.1)) then
		if(anscent.eq.ansbest) then
		rabest=racent
		decbest=deccent
		angbest=angle
		end if
	call dr2tf(3,racent,sign,IDMSt)
	call dr2af(2,deccent,sign,IDMSa)
	write(6,5801) idmst,sign,idmsa,angle,anscent
	write(8,5801) idmst,sign,idmsa,angle,anscent
5801	format(1x,'Original center=',i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x
     .a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,' Ag=',f8.3,
     .' summed weights=',f10.2)
	racent=rabest
	deccent=decbest
	angle=angbest
	call dr2tf(3,racent,sign,IDMSt)
	call dr2af(2,deccent,sign,IDMSa)
	write(6,5802) idmst,sign,idmsa,angle,ansbest
	write(8,5802) idmst,sign,idmsa,angle,ansbest
5802	format(1x,'Revised  center=',i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x
     .a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,' Ag=',f8.3,
     .' summed weights=',f10.2)
	itweak=1
	iaweak=1
	tweakrads=0.0d0
	deltang=0.0d0
	go to 5000
	end if
C
C We are done, except for writing things out!
600	continue
	if(nobjs.gt.0) then
	do 601 i=1,nobjs
	if(objstatus(i).gt.-1) then
C       write(9,1501) i,itrackobjs(i),
C    .                 ixobjs(i),iyobjs(i),ifibs(objstatus(i))
1501	format(1x,'Object',2i5,' at',2i9,' assigned to fiber',i5)
C write exec file:
Cwrite(4,2001) itrackobjs(i),ifibs(objstatus(i))
2001	format('assign',2i9)
	end if
601	continue
	end if
C
	
C
C Start writing things to the output file. Let's start with a warning, in
C case the number of required fops and skys were not met.
        if(nf.lt.2) write(2,5700)
	if(nf.lt.nfopsrequired) write(2,5701)
        if(ns.lt.nskysrequired) write(2,5702)
5700	format('##WARNING! INVALID NUMBER (<2) OF FOPS ASSIGNED!')
5701    format('##WARNING! NUMBER OF FOPS REQUESTED COULD',
     &  ' NOT BE ASSIGNED!')
5702    format('##WARNING! NUMBER OF SKYS REQUESTED COULD',
     &  ' NOT BE ASSIGNED!')
C Now write out the the center & angle
	call dr2tf(3,racent,sign,IDMSt)
	call dr2af(2,deccent,sign,IDMSa)
	write(2,9001)idcent,namecent,IDMSt,sign,IDMSa,angle
9001	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,1x,'C',2x,'Ag=',f8.3)
	if(nacqs.gt.0) then
	itype='A'
	do 701 i=1,nacqs
C write them out, even if we didn't do anything with them
	call dr2tf(3,raacqs(i),sign,IDMSt)
	call dr2af(2,decacqs(i),sign,IDMSa)
Cwrite(3,2002) idacqs(i),nameacqs(i),IDMSt,sign,
C    .                IDMSa
	write(2,9002) idacqs(i),nameacqs(i),IDMSt,sign,
     .                IDMSa,itype
701	continue
	end if
C
C
	if(nfops.gt.0) then
	do 602 i=1,nfops
	call dr2tf(3,rafops(i),sign,IDMSt)
	call dr2af(2,decfops(i),sign,IDMSa)
C write out "leftovers:"
2002	format(3x,i3,4x,a20,5x,i2.2,1x,i2.2,1x,i2.2,'.',i3.3,
     .1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,3x,'Status=',a4)
C figure out the status word
	call itrick(fopstatus(i),ifopsfibs,stat)
C figure out the object type
	itype='F'
Cwrite(3,2002) idfops(i),namefops(i),IDMSt,sign,
C    .                IDMSa,stat
	write(2,9002) idfops(i),namefops(i),IDMSt,sign,
     .                IDMSa,itype,stat
9002	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,1x,a1,2x,'STATUS='
     . ,a4,5x,a1)
	if(fopstatus(i).gt.-1) then
C    	write(9,1502) i,itrackfops(i),
C    .      ixfops(i),iyfops(i),ifopsfibs(fopstatus(i))
1502	format(1x,'Fobs  ',2i5,' at',2i9,' assigned to fiber',i5)
	end if
602	continue
	end if
	if(nskys.gt.0) then
	do 603 i=1,nskys
	call itrick(skystatus(i),ifibs,stat)
	itype='S'
	call dr2tf(3,raskys(i),sign,IDMSt)
	call dr2af(2,decskys(i),sign,IDMSa)
Cwrite(3,2002) idskys(i),nameskys(i),IDMSt,sign,
C    .      IDMSa,stat
	write(2,9002) idskys(i),nameskys(i),IDMSt,sign,
     .      IDMSa,itype,stat
	if(skystatus(i).gt.-1) then
C       write(9,1503) i,itrackskys(i),
C    .                ixskys(i),iyskys(i),ifibs(skystatus(i))
1503	format(1x,'Skys  ',2i5,' at',2i9,' assigned to fiber',i5)
	end if
603	continue
	end if
	if(nobjs.gt.0) then
	do 604 i=1,nobjs
	if(objstatus(i).lt.0) then
	call itrick(objstatus(i),ifibs,stat)
	itype='O'
	call dr2tf(3,raobjs(i),sign,IDMSt)
	call dr2af(2,decobjs(i),sign,IDMSa)
Cwrite(3,2002) idobjs(i),nameobjs(i),IDMSt,sign,
C    . IDMSa,stat
	write(2,9002) idobjs(i),nameobjs(i),IDMSt,sign,
     . IDMSa,itype,stat
C       write(9,1504) i,itrackobjs(i),
C    .ixobjs(i),iyobjs(i),objstatus(i)
1504	format(1x,'Object',2i5,' at',2i9,' was unassigned; status=',i5)
	end if
604	continue
	end if
	if(nobjs.gt.0) then
	do 607 i=1,nobjs
	if(objstatus(i).gt.-1) then
	itype='O'
	call itrick(objstatus(i),ifibs,stat)
	call dr2tf(3,raobjs(i),sign,IDMSt)
	call dr2af(2,decobjs(i),sign,IDMSa)
Cwrite(3,2002) idobjs(i),nameobjs(i),IDMSt,
C    .    sign,IDMSa,stat
	write(2,9002) idobjs(i),nameobjs(i),IDMSt,
     .    sign,IDMSa,itype,stat
	end if
607	continue
	end if
C
	if(nfops.gt.0) then
	do 605 i=1,nfops
Cif(fopstatus(i).lt.0)
C    .  write(9,1505) i,itrackfops(i),
C    .        ixfops(i),iyfops(i),fopstatus(i)
1505	format(1x,'Fobs  ',2i5,' at',2i9,' was unassigned; status=',i5)
605	continue
	end if
	if(nskys.gt.0) then
	do 606 i=1,nskys
Cif(skystatus(i).lt.0)
C    .  write(9,1506) i,itrackskys(i),
C    .   ixskys(i),iyskys(i),skystatus(i)
1506	format(1x,'Sky   ',2i5,' at',2i9,' was unassigned; status=',i5)
606	continue
	end if
C See over what angle this is good for!
C
C
C
C
C
C OK, just about done!
C	write(*,1507)
C	write(8,1507)
C1507	format(1x,'Checking for angle limits: this could take a moment;',
C     .' CNTL-C to skip this step')
C	call anglecheck(ixobjs,iyobjs,nobjs,objstatus,ifibs,
C    .                        ixfops,iyfops,nfops,fopstatus,ifopsfibs,
C     .                        ixskys,iyskys,nskys,skystatus)
C
C
C That's about it, folks!
	stop
	end			
C
C
C
	subroutine anglecheck(ixobjs,iyobjs,nobjs,objstatus,ifibs,
     .                        ixfops,iyfops,nfops,fopstatus,ifopsfibs,
     .                        ixskys,iyskys,nskys,skystatus)
	integer*4 ixobjs(*),iyobjs(*),nobjs,objstatus(*),ifibs(*),
     .            ixfops(*),iyfops(*),nfops,fopstatus(*),ifopsfibs(*),
     .            ixskys(*),iyskys(*),nskys,skystatus(*)
	real*8 angle,PI
	parameter(PI=3.141592653589793d0)
	call rotmode(1)
	isign=1
1	continue
	do 100 i=1,101
	degrees=(i-1)*.01*isign
	angle=degrees*PI/180.
	call rotangle(angle)
Cwrite(9,3001) degrees,angle
3001	format(1x,'Trying angle=',f7.2,2x,g12.6)
C Now move each fiber to its new place:
	if(nobjs.gt. 0) then
		do 101 iobjs=1,nobjs
		if(objstatus(iobjs).gt.-1) then
C it's assigned to a fiber ifibs(objstatus(iobjs))
		call setfib(ifibs(objstatus(iobjs)),ixobjs(iobjs),
     .                       iyobjs(iobjs),istat)
Cwrite(9,3002) iobjs,objstatus(iobjs),ifibs(objstatus(iobjs)),
C    .   ixobjs(iobjs),iyobjs(iobjs),istat
3002	format(1x,'Checking object',i5,' status=',i5,' fiber',i5,
     .  ' x,y=',2i9, 'status=',i5)
			if(istat.ne.0)  then
Cwrite(9,1001) degrees,iobjs,ifibs(objstatus(iobjs)),
C    .               ixobjs(iobjs),iyobjs(iobjs)
1001	format(1x,f7.2,'Angle exceeded for object',i5,' fiber',i5,
     . 'at x,y=',2i9)
			go to 2101
			end if
		end if
101		continue
	end if
	if(nfops.gt. 0) then
		do 102 ifops=1,nfops
		if(fopstatus(ifops).gt.-1) then
C it's assigned to a fiber ifopsfibs(fopstatus(ifops))
		call setfib(ifopsfibs(fopstatus(ifops)),ixfops(ifops),
     .                        iyfops(ifops),istat)
		if(istat.ne.0) then
Cwrite(9,1003) degrees,ifops,ifopsfibs(fopstatus(ifops)),
C    .               ixfops(ifops),iyfops(ifops)
1003	format(1x,f7.2,'Angle exceeded for fops    ',i5,' fiber',i5,
     . 'at x,y=',2i9)
	 go to 2101
	end if
		end if
102	continue
	end if
	if(nskys.gt. 0) then
		do 103 iskys=1,nskys
		if(skystatus(iskys).gt.-1) then
C it's assigned to a fiber ifibs(skystatus(iskys))
		call setfib(ifibs(skystatus(iskys)),ixskys(iskys),
     .                        iyskys(iskys),istat)
		if(istat.ne.0)  then
Cwrite(9,1002) degrees,iskys,ifibs(skystatus(iskys)),
C    .                ixskys(iskys),iyskys(iskys)
1002	format(1x,f7.2,'Angle exceeded for sky   ',i5,' fiber',i5,
     . 'at x,y=',2i9)
	
	go to 2101
	end if
		end if
103	continue
	end if
C OK, well all the fibers REACH the new location. But are there any new
C collisions?
	do 107 ib=1,288
	ibut=ib-1
	itest=icheckall(ibut)
	if(itest.eq.1)  then
Cwrite(9,1004)degrees,ibut
1004	format(1x,f7.2,'Angle exceeded for collision button',i5)
	go to 2101
	end if
107	continue
100	continue
2101	if(isign.eq.1) then
	amax=degrees-.01
	if(amax.lt.0.) amax=0.
	isign=-1
	go to 1
	else
	amin=degrees+.01
	if(amin.gt.0.) amin=0.
	write(2,105) amin,amax
105	format('##Assignments good from angles of',f7.2,' to',f7.2)
	write(*,105) amin,amax
	write(8,105) amin,amax
	return
	end if
	end
C
C
C
	subroutine countem(n,status,nassigned,aweight,wts)
	real*8 wts
	real*4 aweight
	integer*4 status(*)
	nassigned=0
	wts=0.
	if(n.gt.0) then
		do 1 i=1,n
		if(status(i).ge.0) then
		nassigned=nassigned+1
		inum=i
		wts=wts+weight(inum,aweight)
		end if
1		continue
	end if
Cwrite(6,124) n, wts
124	format(1x,'Total in countem for', i5,' objs=', g12.6)
	return
	end

	subroutine itrick(istat,istone,stat)
	integer*4 istone(*)
	character*4 stat,conditions(4)
	data conditions/'OK  ','STP1','STP2','EDGE'/
	if(istat.le.-4) stat='EDGE'
	if(istat.eq.-3) stat='STP2'
	if(istat.eq.-2) stat='STP1'
	if(istat.eq.-1) stat='OK  '
	if(istat.gt.-1) write(stat,'(i4)' ) istone(istat)
	return
	end

C 
C 
	subroutine finalverify(nfibs,ifibs,fibstatus,objstatus,nobjs,
     .             skystatus,nskys,
     .             nfopsfibs,ifopsfibs,fopsfibstatus,fopstatus,nfops)
	integer*4 ifibs(*),fibstatus(*),objstatus(*),ifopsfibs(*),
     .            fopsfibstatus(*),fopstatus(*),skystatus(*)
C this subroutine gently removes anything for which there is still
C a collision.  It should never do anything.
	ieliminate=0
        do 2 i=1,nskys
        if(skystatus(i).gt.-1) then
C it's assigned
        ibut=ifibs(skystatus(i))
        itest=icheckall(ibut)
C       write(9,1112) i,skystatus(i),ifibs(skystatus(i)),itest
1112    format(1x,'final verify checking sky',i5,' stat',i5,' fib',i5, 
     .' test=',i5) 
        if(itest.eq.1) then
C       write(9,1102) i,skystatus(i),ifibs(skystatus(i))
1102    format(1x,'Final verification: eliminating sky',i3,2i5)
                call park(ibut,0)
                ieliminate=ieliminate+1
                fibstatus(skystatus(i))=-1
                skystatus(i)=-1
        end if 
        end if
2       continue
	do 1 i=1,nobjs
	if(objstatus(i).gt.-1) then
C it's assigned
	ibut=ifibs(objstatus(i))
1111	format(1x,'final verify checking obj',i5,' stat',i5,' fib',i5,
     .' test=',i5)
	itest=icheckall(ibut)
Cwrite(9,1111) i,objstatus(i),ifibs(objstatus(i)),itest
	if (itest.eq.1) then
Cwrite(9,1101) i,objstatus(i),ifibs(objstatus(i))
1101	format(1x,'Final verification: eliminating obj',i3,2i5)
		call park(ibut,0)
		ieliminate=ieliminate+1
		fibstatus(objstatus(i))=-1
		objstatus(i)=-1
	end if
	end if
1	continue
	do 3 i=1,nfops
	if(fopstatus(i).gt.-1) then
C it's assigned
	ibut=ifopsfibs(fopstatus(i))
	itest=icheckall(ibut)
cwrite(9,1113) i,fopstatus(i),ifopsfibs(fopstatus(i)),itest
1113	format(1x,'final verify checking obj',i5,' stat',i5,' fib',i5,
     .' test=',i5)
	if(itest.eq.1) then
Cwrite(9,1103) i,fopstatus(i),ifopsfibs(fopstatus(i))
1103	format(1x,'Final verificatin: eliminating fop',i3,2i5)
		call park(ibut,0)
		ieliminate=ieliminate+1
		fopsfibstatus(fopstatus(i))=-1
		fopstatus(i)=-1
	end if
	end if
3	continue
Cwrite(9,1105) ieliminate
1105	format(1x,'Final verification: eliminated total of',i5)
	return
	end
C
C
C
C
	subroutine wings(ixobjs,iyobjs,nthings,status,
     .            nfibers,idfibers,fiberstatus,nlimit)
	integer*4 ixobjs(*),iyobjs(*),status(*),idfibers(*),
     .            fiberstatus(*)
	parameter(PI=3.141592653589793)
	if(nthings.lt.1) return
	if(nfibers.lt.1) return
	if(nlimit.lt.1) return
	nsuccess=0
C OK, now for each object
	do 11 i=1,nthings
	if(nsuccess.eq.nlimit) return
C stop if we've already put enough on.

C find the closest fiber for objects which are as yet unassigned
C and which legally can be assigned (e.g., status=-1)
	if(status(i).eq.-1) then
C	write(9,1101) i,ixobjs(i),iyobjs(i),status(i)
1101		format(1x,'FREEBIE: See if we can sneak in',4i9)
		x=ixobjs(i)
		y=iyobjs(i)
		angle=atan2(y,x)
		if(angle.lt.0.) angle=angle+2.*PI
1		infib=-1
Cwrite(9,1205) (j,idfibers(j),fiberstatus(j),j=1,nfibers)
1205	format(1x,'using fibers',3i5)
		call closestlegalfib(angle,idfibers,nfibers,
     .                              infib,fiberstatus,
     .				    ixobjs(i),iyobjs(i),idone)
C        write(9,1202) idone
1202	format(1x,'Attempt to find closest fiber: successful? 1=yes',i5)
		if(idone.eq.1) then
C See if there are any collisions!
C	write(9,1203) infib,idfibers(infib)
1203		format(1x,'Namely:',2i5)
C Now try moving the fiber there:
		call setfib(idfibers(infib),ixobjs(i),iyobjs(i),istat)
		if(istat.ne.0) stop 'ARHGHWHWHWH!H!!!!!!'
C Now check for collisions:
		itest=icheckall(idfibers(infib))
C	write(9,1204) itest
1204	format(1x,'checking for collisions: 0=none, 1=some',i5)
			if(itest) then
			call park(idfibers(infib),0)
			fiberstatus(infib)=-10
C go try to find another legal fiber, one that won't collide!
			go to 1
			end if
C If we are here then there were no collisions and we were successful!
C set fiber status to a rather special number to show that it is assigned
C in a rather special way---fiberstatus(infib)=0 if it was assigned here
		fiberstatus(infib)=0
		status(i)=infib
		nsuccess=nsuccess+1
Cwrite(9,3001) infib,idfibers(infib),i,ixobjs(i),iyobjs(i)
3001	format(1x,'Sucessfully Snuck On a Freebie: fiber',2i5,'
     .   to',3i9)
		end if
C end the idone=1 loop
C now reset fiberstatusback to -1 from -10
	do 12 j=1,nfibers
	if(fiberstatus(j).eq.-10) fiberstatus(j)=-1
12	continue
	end if 
C end the loop that the object was lega.
11	continue
C do a new object
	return
	end
	function icheckall(ibuta)
2	icheckall=0
	do 1 i=1,288
	ibutb=i-1
	if(ibuta.eq.ibutb) go to 1
	itest=interract(ibuta,ibutb)
	if(itest.eq.1) then
	icheckall=1
Cwrite(9,101) ibuta,ibutb
101	format(1x,'Collision found between buttons',2i6)
	call whereis(ibuta,ix1,iy1,ix2,iy2,ix3,iy3)
Cwrite(9,102) ibuta,ix1,iy1,ix2,iy2,ix3,iy3
102	format(1x,'With locations',i5,' x1,y1:',2i6,' x2,y2:,'2i6,
     .  'x3,y3',2i6)
	call whereis(ibutb,ix1,iy1,ix2,iy2,ix3,iy3)
Cwrite(9,102) ibutb,ix1,iy1,ix2,iy2,ix3,iy3
	return
	end if
1	continue
	return
	end
	
C
C
C
	subroutine gerard(x,y,n,status,nrequired,
     .       nfibass,ifibass,fibasstatus,
     .       nfibs,ifibs,fibstatus,hopstatus,
     .       nfopsfibs,ifopsfibs,fopsfibstatus,aweight)
C this routine takes "n" objects with locations x,y, and status
C and continues to assign them to fibers nfibass,ifibass,fibasstatus
C until all nrequired are done.  They are checked against the "normal"
C fibers nfibs,ifibs,fibstatus, and the fopfibers nfopfibs,ifopfibs,fopsfibstatu
C
C
	parameter(PI=3.141592653589793)
	parameter(maxfibers=100)
	integer*4 x(*),y(*),status(*),ifibass(*),fibasstatus(*),
     .            ifibs(*),fibstatus(*),hopstatus(*),
     .            ifopsfibs(*),fopsfibstatus(*)
C How many of the n objects are currently assigned?
	data piviotfudge/0.0000/
C This is the OLD piviotfudge.
C	data piviotfudge/0.01571/
C
C
Cwrite(9,701)
701	format(1x,'We are now trying to force the issue')
Cwrite(9,702) n,nrequired,nfibass,nfibs,nfopsfibs
702	format(1x,'The number of objects to play with is',i5,
     ./1x,'The number required to be assigned are',i5,
     ./1x,'We can use in principle',i5,' fibers, checked against',2i5)
C Step through the stars:
	if(n.lt.1) return
C
311	icount=0
C
C
C count the number of currently assigned objects and compare it to the 
C required number
	do 305 i=1,n
	if(status(i).gt.-1) icount=icount+1
305	continue
Cwrite(9,1105) icount,nrequired
1105	format(1x,'Now have assigned', i5,' out of',i5,' required')
	if(icount.ge.nrequired) return
C stop when we have done enough!
C
C
C
C
C
C
C See what is the most painless thing.
C initialize pain to the max, and set the "best" fiber and object
C numbers to imaginary values:
	pain=maxfibers
C	write(9,705) pain
705	format(1x,'We are starting with this much pain',f8.1)
	ithisobj=-1
	ithisfib=-1
C
C
C
C
C
C
C Now for each possible object:
	do 306 i=1,n
	if(status(i).ne.-1) go to 306
Cwrite(9,711) i
711	format(1x,' We will see the best we can do with obj',i5)
C Skip if illegal  position or already assigned object
C 
C
C
C Check each of the fibers we can assign to it
C
C
C angle stuff:
C (this is a little obscure.  We need to keep track of the angles
C iff we find two fibers which are equally painless to place on an
C object---we then want to assign the one with the smaller angle.
	angleobj=atan2(float(y(i)),float(x(i)))
	if(angleobj.lt.0.) angleobj=angleobj+2.*PI
C initialize olddif to the largest possible angle
	olddif=PI
C
C
C
C
C
C Now step through each fiber, looking for unassigned fibers that
C will reach object i.
C
C
	do 307 j=1,nfibass
C See if this fiber is sacredly assigned.  If so, skip it:
	if(fibasstatus(j).eq.0) go to 307
	ibut=ifibass(j)
C get the actual fiber number
C Now see if we can assign it:
	call checkfib(ibut,x(i),y(i),istat)
	if(istat.ne.0) go to 307
C skip if if we can't reach it.
Cwrite(9,1106) i,j,ibut,ithisobj,ithisfib,pain
1106	format(1x,'checking object',i5,' with fib',2i5, 
     .'best so far are',2i5,' pain=',
     .f6.2)
C figure out the angle, and the angle difference between this
C fiber and the object:
	anglefib=ibut*2.*PI/288.+piviotfudge
	angledif=abs(anglefib-angleobj)
	if(angledif.gt.PI) angledif=2.*PI-angledif
C (take complement of angle if it is greater than PI!)
C initialize the current pain to no pain at all.
	paina=0.
C
C
C OK, now
C save the position of this fiber:
C note that whereis returns the PHYSICAL coordinate ix1,iy1 and so THIS
C fiber can be returned to this place only by setfibphysphil
	call whereis(ibut,ix1,iy1,ix2,iy2,ix3,iy3)
C	write(9,1414) ibut,ix1,iy1,ix2,iy2
1414	format(1x,'Checking the cell struture:',i5,4i9)
C (ix1,iy1) is where we want to return this fiber!
C make the tenative assignment
	call setfib(ibut,x(i),y(i),istat)
	if(istat.ne.0) stop 'Oh shoot!!!!!!'
C if this fiber was previously assigned, to a non-sacred position,
C note the amount of pain involved:
	if(fibasstatus(j).gt.0) paina=paina+weight(fibasstatus(j),aweight)
Cwrite(9,3000) ibut, paina
3000	format(1x,'starting with this fiber with pain level=',i5,f6.2)
C Now see if there are any collisions.
	itest=icheckall(ibut)
	if(itest.eq.1) then
C there's a collision here, not too surprisingly (but not necessarily!)
C let's see how bad it is!
C
C
C
C First check to see if there is a collision with some fiber other than
C the "current" cable.  This will flag as unacceptable collisions with
C the fops fibers, and collisions with "parked" members of the "other"
C fiber cable.
C  Can only ever happen to the fiber on either side the
C way the angle restriction is now, but lets play it safe! It's only
C CPU cycles!
		do 3002 k=1,288
		ibuta=k-1
C See if this fiber is unassigned or not
		do 3003 ik=1,nfibs
		if(ifibs(ik).eq.ibuta) go to 3002
3003		continue
C ***NOT*** a member of the current fiber bundle:
C
		itest=0
		if(ibuta.ne.ibut)itest=interract(ibuta,ibut)
		if(itest.gt.0) then
		paina=2.*maxfibers
C	write(9,3004) ibut,ibuta
3004	format(1x,'Collision with the OTHER fiber bundle/FOPS',2i5)
		go to 2001
		end if
3002		continue
C OK, there are no collisions with "the other" cable.
C
C Now lets see how bad the collisions are with our own fiber cable
C
C
C
		do 3005 k=1,nfibs
		ibuta=ifibs(k)
		if(ibut.eq.ibuta) go to 3005
C we have already dealt with the situation that the fiber is itself and
C may be assigned.
		itest=interract(ibut,ibuta)
		if(itest.eq.1)  then
C is the fiber ibut is touching (ibuta) assigned to a normal object?
			if(fibstatus(k).lt.1) then
			paina=2.*maxfibers
			go to 2001
C we dont want it!
			end if
C if we UNASSIGN this fiber would it actually do us any good?
C see where the fiber is now...
C note that whereis returns "physical coordinates" ixa1,iya1,
C and so setfibphysphil will be needed to return it to these positions
C (this matters iff concentricity or rotation is on and the fiber is
C parked)
		call whereis(ibuta,ixa1,iya1,ixa2,iya2,ixa3,iya3)
C	write(9,1414) ibuta,ixa1,iya1,ixa2,iya2
C stow the fiber
		call park(ibuta,0)
C		call whereis(ibuta,ixz,iyz,ixy,iyy,ixx,iyx)
C		write(9,1414) ibuta,ixz,iyz,ixy,iyy,ixx,iyx
C now see if there still is a collision
		itest=interract(ibut,ibuta)
C now put it back
C note the use of PHYSICAL coordinates here....
		call setfibphysphil(ibuta,ixa1,iya1,istat)
		if(istat.ne.0) stop 'self-consistency failure test'
		if(itest.eq.1) then
C it's not going to help!
		paina=2.*maxfibers
		go to 2001
		end if
		paina=paina+weight(fibstatus(k),aweight)
		end if
3005		continue
C whew! we have now accumilated a lot of pain.
	end if
C end the conidtional that began with there being some collision
	
2001	continue
C reset fiber ibut to original position
	call setfibphysphil(ibut,ix1,iy1,istat)
Cwrite(9,1415) ibut,ix1,iy1
1415	format(1x,'About to move button',i5,' to',2i9)
	if(istat.ne.0) stop 'what gives? here?'
C skip if paina is eq to pain and the angle is worse
	if((paina.eq.pain).and.(angledif.gt.olddif)) go to 307
C otherwise, if the pain is less than take this one.
	if(paina.le.pain) then
	olddif=angledif
	pain=paina
	ithisobj=i
	ithisfib=j
	end if
307	continue
Cwrite(9,715) i,paina,pain,ithisobj,ithisfib
715	format(1x,'FINALLY: for this object',i5,' the pain level is',
     ./3x,f6.2,' ; the best pain so far is',f6.2,' corrend to obj,fib',
     .2i5)
306	continue
C
C
C
	if(ithisfib.gt.-1) then
Cwrite(9,716) ithisfib,ifibass(ithisfib),ithisobj,
C    .x(ithisobj),y(ithisobj),pain
716	format(1x,'FORCED issue: fib=',2i5,'on',3i7,'pain=',f6.2)
C we were sucessful in finding something that will work
C OK, now the least damaging thing will happen if we assign fops
C thisfib to object ithisobj
		ibut=ifibass(ithisfib)

C So do it!
		call setfib(ibut,x(ithisobj),y(ithisobj),istat)
	if(istat.ne.0) stop 'required assignement failed---why????'
C
C
C See if this fiber was one that was already assigned.
		do 4001 k=1,nfibs
		if(ifibs(k).eq.ibut) then
		if(fibstatus(k).gt.0) then
C fiber is previously assigned and we need to unassign it.
			hopstatus(fibstatus(k))=-1
C unassign the object
			fibstatus(k)=-1
	end if
	if(fibstatus(k).eq.0) then
C something is dreadfully wrong!
	stop 'problem in gerard; fibstatus=0 attempted assignment'
	end if
			go to 4002
		end if
4001		continue
4002		continue
C assign the fiber
		status(ithisobj)=ithisfib
		fibasstatus(ithisfib)=0
C now find the object fibers that conflict with this and undo the assignements
		do 310 k=1,nfibs
		if(fibstatus(k).gt.0) then
			ibuta=ifibs(k)
			itest=0
			if(ibut.ne.ibuta) itest=interract(ibuta,ibut)
			if (itest) then
			call park(ifibs(k),0)
C		write(9,1110) ithisobj,ithisfib,
C    .			ifibass(ithisfib),
C    .    		 k,ifibs(k),fibstatus(k),
C    .			hopstatus(fibstatus(k))
1110	format(1x,'To accomidate obj',i3,' fib',2i5,' we are unassigning',
     ./1x,' normal fib',2i5,' from object',2i5)
			hopstatus(fibstatus(k))=-1
			fibstatus(k)=-1
			end if
		end if
310	continue
	go to 311
	else
Cwrite(9,1112)
1112	format(1x,'We are unable to assign the required number of objs')
	end if
C we now either do or do not have nfopsrequired fibers assigned!
C do the skies!
	return
	end

C This fix introdcued 10/2/95 to provide a little "looseness" in setfibphys
	subroutine setfibphysphil(ibut,ix,iy,istat)	
	call setfibphys(ibut,ix,iy,istat)
	if(istat.eq.0) return
	do 1 i=1,3
	do 1 j=1,3
	ixt=ix+i-2
	iyt=iy+j-2
	call setfibphys(ibut,ixt,iyt,istat)
	if (istat.eq.0) go to 2
1	continue
2	continue
	return
	end
C
	function weight(iobj,aweight)
C	aweight=1.0 is "strong" weighting: 10 objects of weight 10 is
C                                      worth 1 object of weight 1
C              =0.5 is "weak" weighting: 3.333 objects of weight 10 is
C                 worth 1 object of weight 1
C	       =0.0 is "no" weighting: all objects are equal
	weight=(float(iobj))**(-aweight)
	return
	end


	subroutine fake(ixobjs,iyobjs,nobjs)
	integer*4 ixobjs(*),iyobjs(*)
	dimension r(4),theta(4)
	data r/0.5,0.5,0.5,0.5/
	data theta/4.5,8.1,11.7,15.3/
	data parkcircle/76200./
C change from 86300
	nobjs=4
	do 1 i=1,nobjs
	ixobjs(i)=parkcircle*r(i)*cosd(theta(i))+0.5
	iyobjs(i)=parkcircle*r(i)*sind(theta(i))+0.5
1	continue
Cwrite(9,2) (i,ixobjs(i),iyobjs(i),i=1,nobjs)
2	format(1x,'Fake objs:',3i9)
	return
	end
	subroutine fakefops(ixfops,iyfops,nfops,nrequired)
	integer*4 ixfops(*),iyfops(*)
	dimension r(1),theta(1)
	data r/0.5/
	data theta/8./
	data parkcircle/76200./
	nfops=1
	nrequired=1
	do 1 i=1,nfops
	ixfops(i)=parkcircle*r(i)*cosd(theta(i))+0.5
	iyfops(i)=parkcircle*r(i)*sind(theta(i))+0.5
1	continue
Cwrite(9,2) (i,ixfops(i),iyfops(i),i=1,nfops)
2	format(1x,'Fkae fops:',3i9)
	return
	end

	subroutine fakeskys(ixskys,iyskys,nskys,nrequired)
	integer*4 ixskys(*),iyskys(*)
	dimension r(2),theta(2)
	data r/0.1,0.1/
	data theta/180.,100./
	data parkcircle/76200./
	nskys=2
	nrequired=1
	do 1 i=1,nskys
	ixskys(i)=parkcircle*r(i)*cosd(theta(i))+0.5
	iyskys(i)=parkcircle*r(i)*sind(theta(i))+0.5
1	continue
Cwrite(9,2) (i,ixskys(i),iyskys(i),i=1,nskys)
2	format(1x,'Fake skys:',3i9)
	return
	end

	
	
C
C
C
C Preliminary assignment
C
C
C Based closely on Tod Lauer's program as adapted by Taft
C
C
	subroutine tod(xhop,yhop,ihop,hopstatus,nfibs,ifib,fibstatus)
	parameter(maxfibers=100)
	parameter(PI=3.141592653589793)
	integer*4 xhop(*),yhop(*),hopstatus(*),fibstatus(*),ifib(*)
	dimension rad(maxfibers),angle(maxfibers),sbuf(maxfibers)
	dimension key1(maxfibers),key2(maxfibers)
	logical changed
	data parkcircle/76200./
C
C
C
C
C tell it that this fiber is not assigned to an object yet
	do 1 i=1,nfibs
	fibstatus(i)=-1
1	continue
C
C
C find the radius and angle for each object
	do 2 i=1,ihop
	x=xhop(i)
	y=yhop(i)
	rad(i)=1.-sqrt(x**2+y**2)/parkcircle
	angle(i)=atan2(y,x)
	if(angle(i).lt.0.) angle(i)=angle(i)+2.0*PI
	tangle=angle(i)*360./2./PI
Cwrite(9,1101) i,x,y,rad(i),tangle
1101	format(1x,'i,x,y,rad,angle=',i3,2f10.0,2f10.3)
2	continue
1701	format(1x,'returned from sort',i5)
	call sort(ihop,rad,key1)
Cwrite(9,1701) ihop
C	write(9,1102) ihop,(rad(key1(i)),key1(i),i=1,ihop)
1102	format(1x,'ihop,rad,key',i5,/1x,(f10.3,i5))
	call break(rad,key1,ihop,ivalue)
C	write(9,1103) ivalue
1103	format(1x,'ivalue=',i5)
	if(ivalue.gt.0) then
C Begin the loop in which we do the preliminary assignment starting from
C the edge of the plate inwards.  Do this only for the objects that
C are in the outer 3/4rds of the plate, though.  ivalue marks the
C transition
C
	do 110 i=1,ivalue
	keyi=key1(i)
C find the closest free fiber
C jeepers, did Tod ever make heavy going of this in his source code!
C it seems like it becomes EASIER if you do the general case.
C
C
C
C find the unused, legal fiber that is closest in angle 
C	write(9,1106) angle(keyi)
1106	format(1x,'about to go looking for the closestfib for',f8.3)
	call closestlegalfib(angle(keyi),ifib,nfibs,infib,fibstatus,
     .        xhop(keyi),yhop(keyi),idone)
		tangle=360.*angle(keyi)/2./PI
	if(idone.eq.1) then
Cwrite(9,1107) tangle,infib,ifib(infib)
1107	format(1x,'Successfully found a legal fiber for angle',f10.3,
     .': namely,',2i5)
	else
Cwrite(9,1207) tangle
1207	format(1x,'Could not find a legal unused fiber for',f10.2)
	end if
C make the assignment, and flag the fiber as used.
C after up-dating data base IF the move is legal
	if (idone.eq.1) then
C Sucessfully found a legal, unused fiber!
	   call setfib(ifib(infib),xhop(keyi),yhop(keyi),ierr)
	   if(ierr.ne.0) stop 'HELlpppspppp!'
C something serioius is wrong here
		fibstatus(infib)=keyi
		hopstatus(keyi)=infib
Cwrite(9,1108) ifib(infib),xhop(keyi),yhop(keyi),keyi
1108	format(1x, 'set fiber',i5,' to',2i8,' obj=',i5)
C Check to see if the new fiber touches its nearest neighbors
		ie=infib+1
		if(ie.gt.nfibs) ie=1
		is=infib-1
		if(is.lt.1) is=nfibs
		if (fibstatus(ie).gt.-1) call baitandswitch(ie,infib,
     .         xhop,yhop,ifib,fibstatus,hopstatus)
		if (fibstatus(is).gt.-1) call baitandswitch(is,infib,
     .         xhop,yhop,ifib,fibstatus,hopstatus)
	end if
110	continue
C 
C
Cwrite(9,1201) (i,key1(i),xhop(key1(i)),yhop(key1(i)),
C    .  hopstatus(key1(i)),i=1,ivalue)
1201	format(1x,'Object',2i5,' located at',2i7,' is assigned to',i5)
C
	end if
C end the loop in which we have done the assignment by distance
C
C
C
C We have now assigned most of the fibers.  Assign the remaining ones
C using the original ordering in the hopper (assumed to be the priority).
C first re-sort:
C
C
Cwrite(9,1211) ihop,ivalue
1211	format(1x,'ihop,ivalue=',2i5)
	if(ihop.gt.ivalue) then
		do 120 i=ivalue+1,ihop
		sbuf(i-ivalue)=key1(i)
C	write(9,1203) i,sbuf(i-ivalue)
1203	format(1x,'stuffing sbuf:',i3,f6.1)
120		continue
		call sort(ihop-ivalue,sbuf,key2)
Cwrite(9,2005)
2005	format(1x,'returned from sort')
C OK, now do the assignements, using the original priorities
	do 121 i=ivalue+1,ihop
	lookup=key2(i-ivalue)+ivalue
	keyi=key1(lookup)
Cwrite(9,1204) i,lookup,keyi,xhop(keyi),yhop(keyi),rad(keyi)
1204	format(1x,'remaining obj resorted:',5i7,' rad=',f8.2)
C according to Tod, that's our boy.
	call closestlegalfib(angle(keyi),ifib,nfibs,infib,fibstatus,
     .        xhop(keyi),yhop(keyi),idone)
	if(idone.eq.1) then
	tangle=360.*angle(keyi)/2./PI
Cwrite(9,1209) keyi,xhop(keyi),yhop(keyi),
C    . tangle,infib,ifib(infib)
1209	format(1x,'Success assign:',i5,' at',2i6,' theta=',f8.3,' 
     . namely fib',2i5)
	else
Cwrite(9,1208) keyi,xhop(keyi),yhop(keyi),
C    . angle(keyi)
1208	format(1x,'failed to assign fiber for object',3i10,f10.3)
	end if
	if(idone.eq.1) then
C we were sucessful
	call setfib(ifib(infib),xhop(keyi),yhop(keyi),ierr)
	if(ierr.ne.0) stop 'ARRGHHHH!!!!!!!'
		fibstatus(infib)=keyi
		hopstatus(keyi)=infib
C Check to see if the new fiber touches its nearest neighbors
		ie=infib+1
		if(ie.gt.nfibs) ie=1
		is=infib-1
		if(is.lt.1) is=nfibs
		if (fibstatus(ie).gt.-1) call baitandswitch(ie,infib,
     .         xhop,yhop,ifib,fibstatus,hopstatus)
		if (fibstatus(is).gt.-1) call baitandswitch(is,infib,
     .         xhop,yhop,ifib,fibstatus,hopstatus)
	end if
121	continue
	do 1301 i=1,ihop
	tangle=180.*angle(i)/PI
Cwrite(9,1302) i,xhop(i),yhop(i),hopstatus(i),tangle,rad(i)
1302	format(1x,'Prelim sum: obj,x,y,fib,ang,rad=',
     .i5,2i8,i3,2f10.2)
1301	continue
	end if
C end the loop due to ivalue not being equal to ihop.
C

C
C
C
C
C
C
C
C at this point, all objects have been kind of assigned.
C It MIGHT be a good time to re-examine our priorities.....
C
C
C but then, maybe not.
C Let's go looking for collisions.
C
C
C
C
	ihalf=nfibs/2+1
C make this always be at least half of the number of fibers
C
C Now work our way around looking to see if any two fibers touch,
C and if so, attempt to resolve conflicts by swaping the assignments,
C if that would be legal.  If multiple conflicts occur, resolve the highest
C number of diffeerences first.  Repeat this loop as often as you care
C to....
C
C
	do 170 ichl=1,5	
	if(ichl.eq.1) maxch=16
	if(ichl.gt.1) maxch=4
	if(ichl.eq.5) maxch=16
	do 151 nch=1,maxch
	changed=.false.
	do 150 i=1,nfibs
		i1=i
		je=i
		if(fibstatus(je).eq.-1) go to 150
		do 152 j=i+1,i+ihalf
			if(j.gt.nfibs) then
				iefib=j-nfibs
			else
				iefib=j
			end if
			if(fibstatus(iefib).eq.-1) go to 152
C See if the fibers touch, and if so, if switching them would
C in fact be legal:

		itest=iswitchtest(i1,iefib,xhop,yhop,ifib,fibstatus)
C remember that i1 is i; we are just doing it this way in case FORTRAN
C does something funky with array indicies
			if(itest) je=iefib
1401	format(1x, 'In checking fiber',i5,' vs',i5,' we find',
     .i5,' which makes je, iefib=',2i5)
C	write(9,1401) i1,iefib,itest,je,iefib
C
C
152		continue
C
C
C
C
C If we can legally make a swap, do so.
		if(je.ne.i1)  then
		call baitandswitch(i1,je,xhop,yhop,ifib,
     .fibstatus,hopstatus)
		changed=.true.
C if we are calling baitswitch, we are pretty confident that the switch
C will work; if it doesn't, something dreadfully is wrong
C	write(9,1402) ie,je
1402	format(1x,'were switching fibers, and marking changed=true',2i6)
		end if
150	continue
	if(.not.changed) go to 153
C exit the loop if nothing was changed
151 	continue
170	continue
153	continue
C this is a change from how Tod was doing it alright
	do 1501 i=1,ihop
Cwrite(9,1502) i,xhop(i),yhop(i),hopstatus(i),fibstatus(hopstatus(i))
1502	format(1x,'Summary after middle pass:',i4,2i8,i4,i4)
1501	continue
	do 1503 i=1,nfibs
	if(fibstatus(i).eq.-1) go to 1503
Cwrite(9,1504) i,ifib(i),fibstatus(i)
1504	format(1x,'These fibers are assigned:',2i5,' to', i5)
1503	continue
C Now check all the fibers with their neighbors in order
	do 141 nch=1,20
Cwrite(9,1601) nch
1601	format(1x,'Entering final pass:',i5)
	changed=.false.
	do 140 i=1,nfibs
C check to see if the new fiber is touching its nearest neirhbor.
C if so, switch fibers
	i1=i
	if(fibstatus(i1).eq.-1) go to 140
C find the next fiber in use to larger numbers:
	ie=i1
144	ie=ie+1
	if(ie.gt.nfibs) ie=1
	if(fibstatus(ie).eq.-1) go to 144
C See if the fibers touch, and if so, if switching them would
C in fact be legal:
	itest=iswitchtest(i1,ie,xhop,yhop,ifib,fibstatus)
Cwrite(9,1602) i1,ie,itest
1602	format(1x,'Checking if fiber',i5,' touches',i5,' and does it?'
     .,i5)
	if(itest) then
C we can (and should) switch them:
	call baitandswitch(i1,ie,xhop,yhop,ifib,fibstatus,hopstatus)
	changed=.true.
Cwrite(9,1603) i1,ifib(i1),fibstatus(i1),i2,ifib(i2),
C    .fibstatus(i2)
1603	format(1x,'Switching fibers so that',2i5,' =',i5,' and',2i5,
     . ' =',i5)
	end if
140	continue
	if(.not.changed) go to 143
C don't repeat the loop if you didn't change anything
141	continue
143	continue
C
C
C
C
C I'm so confused!  Surely we haven't proven anything yet!
C let us now run through each fiber in turn, see if there
C are any unresolved collisions, and if so,
C remove the lowest priority object.
C We do have to be a little systematic about this
C
C
C
C
C Verifification and removal via priority
	call countem(ihop,hopstatus,nassigned,aweight,wts)
Cwrite(9,5001) nassigned,ihop
5001	format(1x,'Initially we assigned',i5,' objects out of',i5)
C
C
C BIG CHANGE***************************************************
C See if we can eliminate a collision by a complete reassignment:
Cwrite(9,1888)
1888	format('Calling refit')
	call refit(ihop,hopstatus,xhop,yhop,nfibs,fibstatus,ifib)
Cwrite(9,1889)
1889	format('Return from refit')
C 
	call verify(ihop,hopstatus,nfibs,fibstatus,ifib)
	call countem(ihop,hopstatus,nassigned,aweight,wts)
Cwrite(9,5002) nassigned,ihop
5002	format(1x,'After verify pass had',i5,' objects out of',i5)
	do 1801 i=1,ihop
Cwrite(9,1802) i,xhop(i),yhop(i),hopstatus(i),fibstatus(hopstatus(i))
1802	format(1x,'Summary after final pass:',i4,2i8,i4,i4)
1801	continue
	do 1803 i=1,nfibs
	if(fibstatus(i).eq.-1) go to 1803
Cwrite(9,1804) i,ifib(i),fibstatus(i)
1804	format(1x,'These fibers are assigned:',2i5,' to', i5)
1803	continue
	return
	end
C
C
C
C	
	subroutine refit(ihop,hopstatus,xhop,yhop,nfibs,fibstatus,ifib)
	integer*4 hopstatus(*),xhop(*),yhop(*),fibstatus(*),ifib(*)
	if(ihop.eq.1) return
	do 1 i=1,ihop
	if (hopstatus(i).eq.-1) go to 1
	if1=hopstatus(i)
	ibut=ifib(if1)
C skip if no legal fiber assigned
C See if there are any collisions:
Cwrite(9,3) ibut
3	format(1x,'Checking button',i5,' for collisions')
	itest=icheckall(ibut)
	if(itest.eq.1) then
Cwrite(9,11) if1,ibut
11	format('Object',i4,' assigned to fiber',i4, 'has a collision. ',
     &   ' Can we fix?')
C If there's a problem see if we can in fact deal with it simply:
C Can we reassign this object
C to any currently unassigned fiber that would result in there being
C no new collisions?
C find the first unused fiber that can reach it:
C first: where is it?
	ix=xhop(i)
	iy=yhop(i)
	do 2 j=1,nfibs
	if (fibstatus(j).eq.-1) then
C fiber is unassigned
	call checkfib(ifib(j),ix,iy,istat)
	if(istat.eq.0) then
Cwrite(9,12) ifib(j)
12	format(1x,'we COULD get unassigned fiber',i5,' on it instead.')
C YES! We COULD get a different fiber on it!  But would it do
C any good?
	call setfib(ifib(j),ix,iy,istat)
	if(istat.ne.0) stop 'HELlllpppspwpwks'
C But are there ANY collisions with it?
C
C
C Park the original fiber
	call park(ibut,0)
	itest=icheckall(ifib(j))
	if(itest.eq.0) then
C It worked
Cwrite(9,13) 
13	format('YES! Making switch official')
C Note that the original is parked.
	itemp=fibstatus(if1)
	fibstatus(if1)=-1
C Note that the new fiber has been assigned to this object
	fibstatus(j)=itemp
C fiber status has been updated.  Now update object status:
	hopstatus(i)=j
	go to 1
	end if
Cwrite(9,14)
14	format('NO! Moving fibers back!')
	call park(ifib(j),0)
	call setfib(ibut,ix,iy,istat)
	if(istat.ne.0) stop 'Oh, no!  HELP!!!!!'
	end if
	end if
2	continue
Cwrite(9,15) 
15	format('There are no unassigned legal fibers for this')
	end if
1	continue
Cwrite(9,16)
16	format(1x,'We are done with trying to refit')
	return
	end
	subroutine verify(ihop,hopstatus,nfibs,fibstatus,ifib)
	integer*4 hopstatus(*),fibstatus(*),ifib(*)
	ieliminate=0
	if(ihop.eq.1) return
	do 1 i=1,ihop
	i1=i
	if(hopstatus(i1).eq.-1) go to 1
	if1=hopstatus(i1)
	ibut=ifib(if1)
C skip if if there is no fiber assigned
C first check to see if there are ANY collisions
	itest=icheckall(ibut)
	if(itest.eq.1) then
Cwrite(9,1101) i,if1,ibut
1101	format(1x,'Found a collision on verification pass:',3i5)
C there is some collision!
C first see if it is with "the other" fiber cable.  
C  If so, eliminate the assignment.
		do 12 k=1,288
		ibuta=k-1
		do 13 l=1,nfibs
		if(ifib(l).eq.ibuta) go to 12
13		continue
		itest1=0
		if(ibuta.ne.ibut)itest1=interract(ibuta,ibut)
		if(itest1.eq.1) then
C we have a collision with "the other" fiber cable.
Cwrite(9,1102) ibuta
1102	format(1x,'Collision was with The Other cable, fiber=',i5)
			call park(ibut,0)
			hopstatus(i1)=-1
			fibstatus(if1)=-1
			ieliminate=ieliminate+1
			go to 1
		end if
12	continue
	end if
1	continue
C OK, we have now stepped through all the fibers checking to see
C if there was a collision with the "other" fiber cable.  Now what
C about the  normal, assigned cable?  Here we will eliminate the
C lowest priority items:
	do 11 i=1,ihop-1
	i1=i
	if(hopstatus(i1).eq.-1) go to 11
	if1=hopstatus(i1)
C now check for collisions with normal, assigned cables
	do 2 j=i+1,ihop
	i2=j
	if(hopstatus(i2).eq.-1) go to 2
C OK, now see if these guys interact:
	if2=hopstatus(i2)
	itest=interract(ifib(if1),ifib(if2))
	if(itest) then
C simply eliminate whichever of these is lower priority
C which is always i2.
Cwrite(9,1103) j,if2,ifib(if2)
C103	format(1x,'Collision was with normal fiber',3i5)
	fibstatus(if2)=-1
	hopstatus(i2)=-1
	call park(ifib(if2),0)
	ieliminate=ieliminate+1
	end if
2	continue
11	continue
C OK, there is---sadly---still a possibility.  Perhaps an assignment
C coincides with a parked fiber of the SAME cable.  It would have been
C nice to have eliminated these first, but I wasn't smart enough to
C do that.  But let's check them now.
	do 15 i=1,ihop
	if(hopstatus(i).eq.-1) go to 15
	if1=hopstatus(i)
	ibut=ifib(if1)
	itest=icheckall(ibut)
	if(itest.ne.0) then
C eliminate!
		ieliminate=ieliminate+1
C fixed oct 97
		call park(ibut,0)
		hopstatus(i)=-1
		fibstatus(if1)=-1
Cwrite(9,1105)
1105	format(1x,'circumstantial evidence: parked collision')
	end if
15	continue
Cwrite(9,1200) ieliminate
1200	format(1x,'On first verificiation eliminated',i5, 'objects')
	return
	end
C
C			
C		
C
C
C   
C
C
	function iswitchtest(i1,i2,xhop,yhop,ifib,fibstatus)
C this routine tests to see if two fibers ifib(i1) and ifib(i2)
C are touching and if we were to switch them, if that would be legal
C but it differs from baitandswitch in that it doesn't actually switch 'em.
	integer*4 ifib(*),fibstatus(*),xhop(*),yhop(*)
	integer*4 xi1,yi1,xi2,yi2
C see if fiber i1 and i2 touch; if so, see if the switch is legal:
C
C
	iswitchtest=0
C means that they are not touching and/or switching would not be legal or
C would not do any good
C an iswitch=1 implies that switching is a Good Thing
C	write(9,1101) i1,ifib(i1),fibstatus(i1),
C     .xhop(fibstatus(i1)),yhop(fibstatus(i1)),
C    .                i2,ifib(i2),fibstatus(i2),
C    .xhop(fibstatus(i2)),yhop(fibstatus(i2))
1101	format(1x,'Check to see if fiber',2i5,' assigned to', i5, 
     .'at ',2i8,' crosses',
     ./22x,2i5,' assigned to', i5, ' at ',2i8)
	itest=interract(ifib(i1),ifib(i2))
C	if(itest.eq.0) write(9,1102)
1102	format(1x,'They dont.')
	if(itest) then
C		write(9,1103)
1103		format(1x,'They do! Would a switch be legal?')
C switch them
C IF both would be legal
		xi1=xhop(fibstatus(i1))
		yi1=yhop(fibstatus(i1))
		xi2=xhop(fibstatus(i2))
		yi2=yhop(fibstatus(i2))
		call checkfib(ifib(i2),xi1,yi1,istati1)
C	 	write(9,1104) i2,ifib(i2),xi1,yi1,istati1
1104	format(1x,'Attempting to move fib',2i5,' to ',2i8,' returns',i5)
		call checkfib(ifib(i1),xi2,yi2,istati2)
C		write(9,1104) i1,ifib(i1),xi2,yi2,istati2
		if((istati1.eq.0).and.(istati2.eq.0)) then
C OK, it would be legal---but would it do any good?
C We can only answer this by temporaryily changing the database---
C sorry about that!
1105	format(1x,'But would it do any good? 1=no; 0=yes',i5)
		call setfib(ifib(i2),xi1,yi1,istati2)
		call setfib(ifib(i1),xi2,yi2,istati1)
		if((istati2.ne.0).or.(istati1.ne.0)) stop 'HHEEELPPP!'
		itest=interract(ifib(i1),ifib(i2))
C		write(9,1105) itest
		if (itest.eq.0) iswitchtest=1
C yes, it would help!
C
C in any event, move things back!
		call setfib(ifib(i1),xi1,yi1,istati1)
		call setfib(ifib(i2),xi2,yi2,istati2)
		if((istati1.ne.0).or.(istati2.ne.0)) stop 'HHEEELPPP!'
		end if
	end if
	return
	end
	subroutine baitandswitch(i1,i2,xhop,yhop,ifib,fibstatus,hopstatus)
C this routine tests to see if two fibers ifib(i1) and ifib(i2)
C are touching, and switches them IF the switch would be legal,
C AND if switching them would alleviate the problem!
	integer*4 ifib(*),fibstatus(*),hopstatus(*),xhop(*),yhop(*)
	integer*4 xi1,yi1,xi2,yi2
C see if fiber i1 and i2 touch; if so, switch them if the switch is legal:
C
C
Cwrite(9,1101) i1,ifib(i1),fibstatus(i1),
C    .xhop(fibstatus(i1)),yhop(fibstatus(i1)),
C    .                i2,ifib(i2),fibstatus(i2),
C    .xhop(fibstatus(i2)),yhop(fibstatus(i2))
1101	format(1x,'Check to see if fiber',2i5,' assigned to', i5, 
     .'at ',2i8,' crosses with',
     ./22x,2i5,' assigned to', i5, ' at ',2i8)
	itest=interract(ifib(i1),ifib(i2))
C	if(itest.eq.0) write(9,1102)
1102	format(1x,'They dont.')
	if(itest) then
C	write(9,1103)
1103	format(1x,'They do!  Now lets see if we can switch them legally')
C switch them
C IF both would be legal
		xi1=xhop(fibstatus(i1))
		yi1=yhop(fibstatus(i1))
		xi2=xhop(fibstatus(i2))
		yi2=yhop(fibstatus(i2))
		call checkfib(ifib(i2),xi1,yi1,istati1)
C	write(9,1104) i2,ifib(i2),xi1,yi1,istati1
1104	format(1x,'Attempting to move fib',2i5,' to ',2i8,' returns',i5)
		call checkfib(ifib(i1),xi2,yi2,istati2)
Cwrite(9,1104) i1,ifib(i1),xi2,yi2,istati2
		if((istati1.eq.0).and.(istati2.eq.0)) then
C OK, both are legal, so do it:
		call setfib(ifib(i2),xi1,yi1,istati1)
		call setfib(ifib(i1),xi2,yi2,istati2)
		if((istati1.ne.0).or.(istati2.ne.0)) then
Cwrite(9,9101)i2,ifib(i2),xi1,yi1,istati1,i1,
C    .ifib(i1),xi2,yi2,istati2
9101	format('Terrible prob here: ifibno,but=',2i5,'x,y,stat=',3i8,
     . /29x,2i5,8x,3i8)
	stop 'Blehsh!'
C this actually happened once!
		end if
C complete the move
			tmpfibstat=fibstatus(i1)
			fibstatus(i1)=fibstatus(i2)
			fibstatus(i2)=tmpfibstat
			hopstatus(fibstatus(i1))=i1
			hopstatus(fibstatus(i2))=i2
C		write(9,1105) i1,ifib(i1),fibstatus(i1),
C     .		xhop(fibstatus(i1)),yhop(fibstatus(i1))
C		write(9,1105) i2,ifib(i2),fibstatus(i2),
C    .		xhop(fibstatus(i2)),yhop(fibstatus(i2))
1105	format(1x,'Fib',2i5,' now assigned to obj',i6,' at',2i8)
1106	format(1x,'But did it actually improve anything? 1=no; 0=yes',i5)
		itest=interract(ifib(i1),ifib(i2))
C		write(9,1106) itest
		if(itest) then
C no, it's still as bad.  Move 'em back!
		call setfib(ifib(i1),xi1,yi1,istati1)
		call setfib(ifib(i2),xi2,yi2,istati2)
		if((istati1.ne.0).or.(istati2.ne.0))stop'HELPF!' 
		tmpfibstat=fibstatus(i1)
		fibstatus(i1)=fibstatus(i2)
		fibstatus(i2)=tmpfibstat
		hopstatus(fibstatus(i1))=i1
		hopstatus(fibstatus(i2))=i2
C		write(9,1105) i1,ifib(i1),fibstatus(i1),
C    .		xhop(fibstatus(i1)),yhop(fibstatus(i1))
C		write(9,1105) i2,ifib(i2),fibstatus(i2),
C    .		xhop(fibstatus(i2)),yhop(fibstatus(i2))
		end if
		end if
	end if
	return
	end


C	
C
C
C this routine checks through all nfibs fibers until it finds the
C unused, legal fiber that has the smallest angle difference between
C its location and the object position
	subroutine closestlegalfib(angle,ifib,nfibs,infib,fibstatus,
     .ix,iy,idone)
C find the closest free fiber
	integer*4 ifib(*),fibstatus(*)
	parameter(PI=3.141592653589793)
C the piviotfudge takes into account that the center of the fibers aren't in
C the right place.
C except that now they are
C	data piviotfudge/0.01571/
C 2PI/400
	data piviotfudge/0.00000/
	infib=-1
	idone=0
	small=2.*PI
Cwrite(9,1100) nfibs,angle
1100	format(1x,'about to go looking through',i5, ' fibers for closest',
     ./' to', f10.3, ' in angle')
	if(nfibs.lt.1) return
	do 1 i=1,nfibs
	angfib=ifib(i)*2.*PI/288.+piviotfudge
	dif=abs(angfib-angle)
C if more than half-way around, take the complement.
	if(dif.gt.PI) dif=2.*PI-dif
	if(fibstatus(i).eq.-1) then
	if(dif.lt.small) then
Cwrite(9,1101) i,ifib(i),angfib,dif
1101	format(1x,'this one is better: item',i5,' fiber',i5,' ang=',2f10.3)
C check to see if the move is legal or not
	call checkfib(ifib(i),ix,iy,istat)
Cwrite(9,1102) ifib(i),ix,iy,istat
1102	format(1x,'But is it legal? Fiber',i5,' move to',2i10,' is',i5)
		if(istat.eq.0) then
			idone=1
			small=dif
			infib=i
		end if
	end if
	end if
1	continue
	return
	end
C
C
C
C	
C If all fibers were present and uniformly spaced, this function would
C be trivial and equal to something like 288*angle/(2*PI)
C
C
	subroutine break(rad,key1,ihop,ivalue)
C ...at a distance of 0.77 inwards we can suddenly reach a given object
C with nearly all (44/50) of the fibers we could eventually reach it
C with.  My guess is that this is a good way to break where we are
C doing the assignment via radius, although more tests would clearly
C be Really Good.
	dimension rad(*)
	dimension key1(*)
	data a0/0.77/
	ivalue=0
	do 1 i=1,ihop
Cwrite(9,1101) i,key1(i),rad(key1(i)),a0,ivalue
1101	format(1x,'trying to find ivalue: i,key,rad,a0,ivalue',2i5,
     .f5.3,f5.3,i5)
	if(rad(key1(i)).gt.a0) return
	ivalue=i
1	continue
	return
	end
C
C
C
C This is stolen directly from numerical recipies
C p. 233 INDEXX
C
	subroutine sort(n,arrin,indx)
	dimension arrin(*),indx(*)
	do 11 j=1,n
	indx(j)=j
11	continue
	if(n.eq.1) return
	l=n/2+1
	ir=n
10	continue
		if(l.gt.1) then
			l=l-1
			indxt=indx(l)
			q=arrin(indxt)
		ELSE
			indxt=indx(ir)
			q=arrin(indxt)
			indx(ir)=indx(1)
			ir=ir-1
			if(ir.eq.1) then
				indx(1)=indxt
				return
			endif
		endif
	i=l
	j=l+l
20	if(j.le.ir) then
		if(j.lt.ir) then
			if(arrin(indx(j)).lt.arrin(indx(j+1)))j=j+1
		endif
		if(q.lt.arrin(indx(j)))then
			indx(i)=indx(j)
			i=j
			j=j+j
		else
			j=ir+1
		endif
	go to 20
	endif
	indx(i)=indxt
	go to 10
	end
	

C
C
C
	subroutine initallfibers(icable,nfibs,ifibs,nfopsfibs,ifopsfibs)
	dimension ifibs(*),ifopsfibs(*)
	character*80 inline
C Major rewrite 19 Nov 1999
C Now uses the local "concentricities" file to determine active fibers
C First, make SURE concentricity is on, rotation is off.
	call concenmode(1)
	call rotmode(0)
C zero the counter for number of fops fibers
	nfopsfibs=0
C zero the counter for number of object fibers
	nfibs=0
C Open the concentricities file
	open(unit=55,file='concentricities',status='old',err=991)
C Now find the right place in the file
1000	read(55,2000,end=992) inline
2000	format(a80)
	if(inline(1:5).ne.'#FIB ') go to 1000
C skip over the dummy line
	read(55,2000) inline
C
C
C
C Now read in each fiber and see if its from the cable we're using
C And whether or not the fiber is "active"
C
C
C
	icheck=-1
2001	read(55,2000,end=999) inline
	icheck=icheck+1
	read(inline(1:3),*,err=993) num
C Make sure that the fiber number is valid
	if(num.ne.icheck) go to 993
	if(num.gt.287) go to 993
C Check to make sure the concentricities file hasn't been messed up by a tab
C character or something equally silly
	if((inline(17:17).ne.'B').and.(inline(17:17).ne.'R').and.
     &     (inline(17:17).ne.'Z').and.(inline(17:17).ne.'F').and.
     &     (inline(17:17).ne.'P').and.(inline(17:17).ne.'S')) go to 994
	if((inline(20:20).ne.'A').and.(inline(20:20).ne.'B').and.
     &     (inline(20:20).ne.'L').and.(inline(20:20).ne.'Z')) go to 995
C OK, this line looks OK.  Now see if it is a cable we want.
C Is it a fobs?
	if((inline(17:17).eq.'F').and.(inline(20:20).eq.'A')) then
	nfopsfibs=nfopsfibs+1
	ifopsfibs(nfopsfibs)=num
	go to 2001
	end if
	if((inline(17:17).eq.'R').and.(inline(20:20).eq.'A').and.
     &   (icable.eq.1))	then
	nfibs=nfibs+1
	ifibs(nfibs)=num
	go to 2001
	end if
	if((inline(17:17).eq.'B').and.(inline(20:20).eq.'A').and.
     &   (icable.gt.1)) then
	nfibs=nfibs+1
	ifibs(nfibs)=num
	go to 2001
	end if
	go to 2001
991	stop 'Cannot find concentricities file.  I quit'
992	stop 'Premature end of concentricities file. I quit'
993	stop 'Problem with format of concentricites file. I quit.'
994	stop 'Invalid cable code in concentricities file. I quit.'
995	stop 'Invalid status code in concenticities file. I quit.'
999	continue
C Report number of fibers found:
	write(6,9101) nfopsfibs
9101	format('There are ',i2,' active fops fibers')
	write(6,9102) nfibs
9102	format('There are ',i3,' active fibers in the',$)
	if(icable.eq.1) then
	write(6,9103)
	else
	write(6,9104)
	end if
9103	format(' RED cable')
9104	format(' BLUE cable')
C Now make sure that all buttons are 'parked'
	do 31 i=1,288
	ibut=i-1
	call park(ibut,0)
31	continue
	close(unit=55)
	return
	end
C This computes the best rotator angle for a given declination
	subroutine getangle(deccent,angle)
	implicit real*8 (a-h,o-z)
	parameter(PI=3.141592653589793d0)
	data alat/31.9633/
	ZD=dacosd(dsin(deccent)*dsind(alat)+dcos(deccent)*dcosd(alat))
	decd=deccent*180./PI
C these two lines changed 17 March 1995
	angle=90.-ZD
	if(decd.gt.alat) angle=-90.-ZD
	return
	end
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
Cwrite(9,9205) xpp,ypp
9205	format(10x,'Encoder steps=',2f10.0/)
C
C
	ix(i)=xpp+0.5
	iy(i)=ypp+0.5
	call legal(xpp,ypp,istat)
	status(i)=istat
C  
C status= -4 beyond park radius
C status= -1 OK.  Use of -1 rather than 0 as we will use this
C     		array for assignment status, and it will thus retain this silly
C     		zero-based numbering scheme for the buttons.
100	continue
	return
	end
C
C
C
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
        subroutine legal(xp,yp,istat)
C modified Aug 10, 1997 to check for objects outside of the parkcircle
C modified Aug 12, 1997 by R. Wolff to start with the closest fiber-10.
        implicit real*8 (a-h,o-z)
        parameter(PI=3.141592653589793d0)
        integer*4 ix,iy,itest
C Hey, is this legal or what?
        data parkcircle/76200./
        data fudge/100./
        istat=-1
C no problem!
C First, see if the position is well beyond the parkcircle:
        if((xp.eq.0.).and.(yp.eq.0.)) go to 112
        radius=sqrt(xp**2+yp**2)
        if(radius.gt.parkcircle+fudge) then
        istat=-4
        return
        end if
112     continue
C next, see if some fiber somewhere can reach this position.  If so, then
C it is OK, and we can leave this routine
C
C       Estimate the number of the closest fiber
        tangle = atan2(yp,xp)
        if (tangle .lt. 0.) tangle = tangle + 2.*PI
C       and back off some
        offset = int(288. * tangle/(2.*PI)) -10
        ix=xp
        iy=yp
        do 1 i=1,288
        ibut= offset + i - 1
C       Modulo arithmetic to keep ibut in range 0 - 287
        if (ibut .lt. 0) ibut = ibut + 288
        if (ibut .ge. 288) ibut = ibut - 288
C       write(6,111) ibut,ix,iy
111     format(1x,'ibut=',3i9)
        call checkfib(ibut,ix,iy,itest)
C       write(6,111) ibut
        if(itest.eq.0) return
1       continue
C mark flag to indicate that this position is illegal, but we don't know why
        istat=-5
C if we made it this far, then it can't reach this position for SOME reason.
        radius=sqrt(xp**2+yp**2)
        if(radius.gt.parkcircle) then
                istat=-4
        end if
        if(istat.eq.-5) write(9,11) ibut,ix,iy,itest
11      format('MYSTERIOUS PROBLEM: ibut, ix,iy, itest=',4i8)
        return
        end

C subroutine for getting the ra and dec out of the file
C
C
	subroutine grabradec(racent,deccent,epochcent,idcent,namecent,
     .    rafops,decfops,epochfops,idfops,namefops,Nfops,itrackfops,
     .     maxfops,
     .    raacqs,decacqs,epochacqs,idacqs,nameacqs,Nacqs,itrackacqs,
     .     maxacqs,
     .	  raskys,decskys,epochskys,idskys,nameskys,Nskys,itrackskys,
     .     maxskys,
     .	  raobjs,decobjs,epochobjs,idobjs,nameobjs,Nobjs,itrackobjs,
     .     maxobjs,
     .                 deepoch,ifile)
	implicit real*8 (a-h,o-z)
	real*8 raobjs(*),decobjs(*),epochobjs(*),
     .         raskys(*),decskys(*),epochskys(*),
     .      rafops(*),decfops(*),epochfops(*),
     .      raacqs(*),decacqs(*),epochacqs(*)
	integer*4 idobjs(*),idskys(*),idfops(*),
     .            idacqs(*),itrackfops(*),itrackacqs(*),itrackskys(*),
     .            itrackobjs(*)
        character*20 nameobjs(*),nameskys(*),
     .       nameacqs(*),namefops(*),namecent
	character*20 honiker
	character*80 header
	character*1 sign
	character*1 class
C
C
C
C
C
C initialize counters
	nobjs=0
	nfops=0
	nacqs=0
	nskys=0
	icount=-1
	if(ifile.eq.1) go to 1000
C We deal now with a Nessie file
C
C
C
C
C Skip over header
	rewind(1)
	read(1,'(A)') header
1	read(1,2,end=3) i,honiker,ihr,imn,sec,sign,IDEG,IMIN,ARCSEC
	if((i.eq.0).and.(honiker.eq.'                    ').and.
     .  (ihr.eq.0).and.(imn.eq.0).and.(sec.eq.0.).and.(sign.eq.' ').and.
     .  (ideg.eq.0).and.(imnin.eq.0).and.(arcsec.eq.0.)) go to 1
C skip over blank line
C	write(9,133) i,honiker
C133	format(1x,i3,2x,'name=',a20)
2	format(3x,i3,4x,a20,5x,i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2)
        ra=STOR(ihr,imn,sec)
        dec=ASTOR(SIGN,IDEG,IMIN,ARCSEC)
	icount=icount+1
C count the number of objects, staring with 0 as the first one.
C this is so that 'assign' exec file will work
	if((honiker(1:6).eq.'Center').or.
     .     (honiker(1:6).eq.'CENTER').or.
     .     (honiker(1:6).eq.'center')) then
C		write(9,1234) ihr,imn,sec,sign,ideg,imin,arcsec
C1234	format(1x,'ihr,imn,sec,etc',1x,i2,1x,i2,1x,f8.3,1x,a1,
C    .i2,1x,i2,1x,f6.2)
		racent=ra
		deccent=dec
		epochcent=deepoch
		idcent=i
		namecent=honiker
		ncent=1
	ELSE if((honiker(1:7).eq.'Acquire').or.
     .          (honiker(1:7).eq.'ACQUIRE').or.
     .          (honiker(1:7).eq.'acquire')) then
		nacqs=nacqs+1
		If (nacqs.gt.maxacqs) then
		write(6,601) maxacqs
601		format(1x,'Number of acq stars exceeds',
     .          ' maximum number', i5'. I quit.')
		stop
		end if
		raacqs(nacqs)=ra
		decacqs(nacqs)=dec
		epochacqs(nacqs)=deepoch
		idacqs(nacqs)=i
		nameacqs(nacqs)=honiker
		itrackacqs(nacqs)=icount
	ELSE if((honiker(1:5).eq.'Guide').or.
     .          (honiker(1:5).eq.'GUIDE').or.
     .		(honiker(1:5).eq.'guide')) then
		nfops=nfops+1
		If (nfops.gt.maxfops) then
		write(6,602) maxfops
602		format(1x,'Number of fops stars exceeds',
     .          ' maximum number', i5'. I quit.')
		stop
		end if
		rafops(nfops)=ra
		decfops(nfops)=dec
		epochfops(nfops)=deepoch
		idfops(nfops)=i
		namefops(nfops)=honiker
		itrackfops(nfops)=icount
	ELSE if((honiker(1:3).eq.'Sky').or.
     .		(honiker(1:3).eq.'SKY').or.
     .  	(honiker(1:3).eq.'sky')) then
		nskys=nskys+1
		If (nskys.gt.maxskys) then
		write(6,603) maxskys
603		format(1x,'Number of skys exceeds',
     .          ' maximum number', i5'. I quit.')
		stop
		end if
		raskys(nskys)=ra
		decskys(nskys)=dec
		epochskys(nskys)=deepoch
		idskys(nskys)=i
		nameskys(nskys)=honiker
		itrackskys(nskys)=icount
	ELSE
		nobjs=nobjs+1
		If (nobjs.gt.maxobjs) then
		write(6,604) maxobjs
604		format(1x,'Number of objects exceeds',
     .          ' maximum number', i5'. I quit.')
		stop
		end if
		raobjs(nobjs)=ra
		decobjs(nobjs)=dec
		epochobjs(nobjs)=deepoch
		idobjs(nobjs)=i
		nameobjs(nobjs)=honiker
		itrackobjs(nobjs)=icount
	
	END IF
	go to 1
3	Continue
C finished reading Nessie file
C Check to see that we do have the plate center!
	if(ncent.ne.1) stop 'No plate center in Nessie file!'
C
C OK, now a little clean up:
C	
C	check to see that there is at least one acquisition star.  If
C       not, pretend there is one at the center!
C This code temporarily removed while I think about it.
C	if(nacqs.eq.0) then
C	nacqs=1
C	raacqs(nacqs)=racent
C	decacqs(nacqs)=deccent
C	epochacqs(nacqs)=epochcent
C	idacqs(nacqs)=idcent
C	nameacqs(nacqs)=namecent
C	end if
	return
C
C
C
C
C
C
C Come here if it is a New Fangled Hydra file
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=1003) 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(2,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.  Let's see if we can read the damn thing.
	read(header,1002) i,honiker,ihr,imn,sec,sign,IDEG,IMIN,ARCSEC,
     .  class
	epoch=deepoch
Cwrite(9,1002) i,ihoniker,ihr,imn,sec,sign,IDEG,IMIN,ARCSEC,
C    .  class
1002	format(i4,1x,a20,1x,i2,1x,i2,1x,f6.3,1x,a1,
     .  i2,1x,i2,1x,f5.2,1x,a1)
	if((i.eq.0).and.(honiker.eq.'                    ').and.
     .  (ihr.eq.0).and.(imn.eq.0).and.(sec.eq.0.).and.(sign.eq.' ').and.
     .  (ideg.eq.0).and.(imnin.eq.0).and.(arcsec.eq.0.)) go to 1001
C skip over blank line
        ra=STOR(ihr,imn,sec)
        dec=ASTOR(SIGN,IDEG,IMIN,ARCSEC)
	icount=icount+1
C count the number of objects, staring with 0 as the first one
	if((class.eq.'C').or.(class.eq.'c'))then
		racent=ra
		deccent=dec
		epochcent=epoch
		idcent=i
		namecent=honiker
		ncent=1
	ELSE if((class.eq.'A').or.(class.eq.'a')) then
		nacqs=nacqs+1
		If (nacqs.gt.maxacqs) then
		write(6,601) maxacqs
		stop
		end if
		raacqs(nacqs)=ra
		decacqs(nacqs)=dec
		epochacqs(nacqs)=epoch
		idacqs(nacqs)=i
		nameacqs(nacqs)=honiker
		itrackacqs(nacqs)=icount
	ELSE if ((class.eq.'F') .or.(class.eq.'f'))then
		nfops=nfops+1
		If (nfops.gt.maxfops) then
		write(6,602) maxfops
		stop
		end if
		rafops(nfops)=ra
		decfops(nfops)=dec
		epochfops(nfops)=epoch
		idfops(nfops)=i
		namefops(nfops)=honiker
		itrackfops(nfops)=icount
	ELSE if((class.eq.'S').or.(class.eq.'s')) then
		nskys=nskys+1
		If (nskys.gt.maxskys) then
		write(6,603) maxskys
		stop
		end if
		raskys(nskys)=ra
		decskys(nskys)=dec
		epochskys(nskys)=epoch
		idskys(nskys)=i
		nameskys(nskys)=honiker
		itrackskys(nskys)=icount
	ELSE
		nobjs=nobjs+1
		If (nobjs.gt.maxobjs) then
		write(6,604) maxobjs
		stop
		end if
		raobjs(nobjs)=ra
		decobjs(nobjs)=dec
		epochobjs(nobjs)=epoch
		idobjs(nobjs)=i
		nameobjs(nobjs)=honiker
		itrackobjs(nobjs)=icount
	
	END IF
	go to 1001
1003	Continue
C finished reading the New Fangled Hydra File
C Check to see that we do have the plate center!
	if(ncent.ne.1) stop 'No plate center in Hydra file!'
C
C OK, now a little clean up:
C	
C	check to see that there is at least one acquisition star.  If
C       not, pretend there is one at the center!
C temporarily removed this code while I thought about it.
C	if(nacqs.eq.0) then
C	nacqs=1
C	raacqs(nacqs)=racent
C	decacqs(nacqs)=deccent
C	epochacqs(nacqs)=epochcent
C	idacqs(nacqs)=idcent
C	nameacqs(nacqs)=namecent
C	end if
	return
C
C
	
	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(2,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)
	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
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 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 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
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   


