! @(#)sparlax.prg	14.1.1.1 (ES0-DMD) 09/16/99 09:38:31
!======================================================================
! MIDAS PROCEDURE   parlax.prg        (A. Smette, April 1994)
!
! Compute the parallactic angle and the atmospheric differential refraction
! cf. Filippenko, A.V.: 1982, PASP, 94, 715
!
!  input data: right ascension in hour,min,seconds (seconds can be fractionary)
!              declination     in degree,min,seconds (idem)
!              sidereal time   in hour,min,seconds  (idem)
!           working wavelength in Angstroms (default: 5000)
!         reference wavelength in Angstroms (default: 5000)
!
!! Note: ra and dec should be apparent coordinates...
!
! assumes std conditions for the refraction angle calculation: 
!  Alt = 2km; phi = +/- 30deg; P = 600 mm Hg; T = 7deg C; f = 8 mm Hg
! (not much dependence on f or T)
!
!  output:      altitude in degrees ( 0 at horizon, 90 at zenith)
!
!               azimuth in degrees (south=0,west=90,north=180,east=270)
!
!		parallactic angle in degrees (-180 to 180)
!
!		differential refraction angle in arcseconds (negative values
!             mean that the working wavelength image is at altitude
!             larger than reference wavelength image
!
! 
! Version 1.0, A. Smette     
!         1.1  Corrected bug April 94; A. Smette
!         1.2  Corrected bug April, 9, 94; A. Smette (thanks to Gerry Williger)
!         1.3  Corrected bug April,11, 94; A. Smette (thanks to Gerry Williger)
!
! example: 
! @@ pa 11,04,0 -18,05,0 11,0,0 3500. 6000.
! Altitude                      : 78.62 degrees 
! Azimuth                       : 184.83 degrees 
! Parallactic Angle             : -175.57 degrees 
! Differential Refraction Angle : -0.31 arcseconds 
!
!----------------------------------------------------------------------
!
! Definition of the parameters
!
define/par p1 ? n "alpha: "
define/par p2 ? n "delta: "
define/par p3 ? n "sideral time: "
define/par p4 5000. n "working   wavelength: "
define/par p5 5000. n "reference wavelength: "
!
! some constants
!
define/local alphac/c/1/10 {p1}
define/local deltac/c/1/10 {p2}
define/local stc/c/1/8 {p3}

define/local indx/i/1/1 0
define/local i/i/1/1 0

! phi is the latitude of La Silla observatory
! (Schmidt telescope: -29deg15'25.8"
!
define/local phi/r/1/1     -29.257167
define/local sinphi/r/1/1  -0.488730271
define/local cosphi/r/1/1   0.872434883
!define/local phi/r/1/1    -29.
!define/local sinphi/r/1/1 -0.48480962
!define/local cosphi/r/1/1  0.87461971
!
! standard atmospheric conditions:
!
define/local p/d/1/1 600.d+00
define/local t/d/1/1 7.d+00
define/local f/d/1/1 8.00d+00

!
! formula are for lambda in microns, so convert them:
!
define/local tmp1/d/1/1 {p4}
define/local tmp2/d/1/1 {p5}
define/local lc/d/1/2 0.d+00,0.d+00
lc(1) = tmp1/10000.d+00
lc(2) = tmp2/10000.d+00
!
define/local l/d/1/1 0.d+00
define/local l2/d/1/1 0.d+00
define/local nc/d/1/2 0.d+00,0.d+00
define/local n/d/1/1 0.d+00
define/local r/d/1/1 0.d+00
define/local z/d/1/1 0.d+00

define/local alpha/r/1/3 0.,0.,0.
define/local delta/r/1/3 0.,0.,0.
define/local st/r/1/3 0.,0.,0.

define/local alphadec/r/1/1 0.
define/local deltadec/r/1/1 0.
define/local stdec/r/1/1 0.


define/local ratio/r/1/1 0.
define/local absratio/r/1/1 0.
define/local hamax/r/1/1 0.

define/local ha/r/1/1 0.
define/local haabs/r/1/1 0.
define/local altitude/r/1/1 0.
define/local azimuth/r/1/1 0.
define/local pa/r/1/1 0.


indx = 0
do i = 1 3
  indx = m$index(alphac,",")
  if indx .eq. 0 goto compute0
  indx = indx-1
  alpha({i}) = {alphac(1:{indx})}
  indx = indx+2
  alphac = "{alphac({indx}:10)}"
enddo
compute0:

indx = 0
do i = 1 3
  indx = m$index(deltac,",")
  if indx .eq. 0 goto compute1
  indx = indx-1
  delta({i}) = {deltac(1:{indx})}
  indx = indx+2
  deltac = "{deltac({indx}:10)}"
enddo
compute1:

indx = 0
do i = 1 3
  indx = m$index(stc,",")
  if indx .eq. 0 goto compute2
  indx = indx-1
  st({i}) = {stc(1:{indx})}
  indx = indx+2
  stc = "{stc({indx}:8)}"
enddo
compute2:

alphadec = (alpha(1)+alpha(2)/60.+alpha(3)/3600.)*15.
deltadec = delta(1)+delta(2)/60.+delta(3)/3600.
stdec = (st(1)+st(2)/60.+st(3)/3600.)*15.

ha = stdec-alphadec
haabs = m$abs(ha)

ratio = m$tan(phi)/m$tan(deltadec)
absratio = m$abs(ratio)

if {absratio} .gt. 1. then 
	hamax = 1000.
else
	hamax = m$acos(ratio)
endif

altitude = m$sin(deltadec)*sinphi+m$cos(deltadec)*m$cos(haabs)*cosphi
pa = altitude
altitude = m$asin(altitude)
z = 90.-altitude

!
! azimuth is computed from south to west north and east
! between 0 and 360 degrees
!
! here, we work with ha always positiv => 0 < azimuth < 180
! 
azimuth = m$cos(haabs)*sinphi-m$tan(deltadec)*cosphi
azimuth = m$sin(haabs)/azimuth
azimuth = m$atan(azimuth)
if azimuth .lt. 1.e-12 azimuth = 180.+azimuth

pa = 1.-pa*pa
pa = m$sqrt(pa)
pa = m$sin(haabs)*cosphi/pa
pa = m$asin(pa)

if deltadec .lt. phi pa = (180.-pa)
if haabs .lt. hamax pa = (180.-pa)
if pa .gt. 180. pa = (180.-pa)
if pa .gt. 360 pa = pa-360.
if pa .lt. -0 pa = pa+360.
if pa .lt. -360 pa = pa+360.

do i = 1 2
 l = lc({i})
 l2 = l*l
 goto refraction
 back:
 nc({i}) = n
enddo

r = 206265.d+00*(nc(2)-nc(1))*m$tan(z)

! for objects east of meridian, change the sign... 
!
if ha .lt. 0. then
  pa = -1.*pa
  azimuth = 360.-azimuth
endif

set/format F6.2
write/out "Altitude                      : {altitude} degrees"
write/out "Azimuth                       : {azimuth} degrees"
write/out "Parallactic Angle             : {pa} degrees"
write/out "Differential Refraction Angle : {r} arcseconds"
set/format

return

refraction:

set/format G19.10
n = 64.328d+00
n = n+29498.1d+00/(146.d+00-1.d+00/l2)
n = n+255.4d+00/(41.d+00-1.d+00/l2)
n = n*p*(1.d+00+(1.049d+00-0.0157d+00*t)*1.d-06*p)
n = n/720.883d+00/(1.d+00+0.003661d+00*t)
n = n-(0.0624d+00-0.000680d+00/l2)/(1.d+00+0.003661d+00*t)*f
n = n*1.d-06
!
! n is in fact n-1
! however, we make the difference between two n so n1-n2 = (n1-1)-(n2-1)
! 
goto back