	implicit real*8 (a-h,o-z)
	A=0.224645E-03
	B=-.257325E-06
	raref=0.204691
	decref=0.530353
	alst=5.7
	exp=0.1
	call antirefract(raref,decref,alst,exp,a,b,ran,dcn)
	write(6,12) raref,decref,ran,dcn
12	format('FINAL',4g16.8)
	stop
	end

	subroutine antirefract(raout,decout,alst,exp,a,b,rain,decin)
C Lets just do it the silly way...
  	implicit real*8 (a-h,o-z)
	inum=0
	rain=raout
	decin=decout
	write(8,*) raout,decout,alst,exp,a,b
12	call refractcor(rain,decin,alst,exp,a,b,ratest,dectest,0)
	inum=inum+1
	difra=(ratest-raout)*206265.
	difdec=(dectest-decout)*206265.
	write(8,13) inum,difra,difdec
13	format(i5,2f10.2)
	if((abs(difra).lt.0.01).and.(abs(difdec).lt.0.01)) return
	rain=rain-difra/2./206265.
	decin=decin-difdec/2./206265.
	go to 12
	end
C 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 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

