;
; Written by Enrico Marchetti, ESO, January 2001
; The routine dar.pro is NOT the main routine.
; The main routine is diff_atm_refr.pro
;
; By the way the wavelength range, the airmass range and site atmospheric parameters
; can be easily modified editing directly this file, saving and recompiling it:
;
; IDL>.r diff_atm_refr
;

Pro DAR,Z,DR,Lambda,Lambda0

; The routines computes the Differential Atmospheric Dispersion
; for a given zenithal distance "Z" for different wavelengths "Lambda"
; with respect to a reference wavelength "Lambda0".
;
; The atmospheric parameters can be adjusted to those characterstic
; of the site the computation is made for.
; The parameters listed below refer to the average Paranal conditions.
;
; CALLING SEQUENCE
;             
;       DAR,Z,DR [,Lambda]
;
;
; ARGUMENTS
;
;       Z 
;       The zenithal distance in degrees
;
; OUTPUT
;
;       DR
;       The Differential Atmospheric Dispersion in arcseconds
;       with respect to the reference wavelength Lambda0
;
;       Lambda
;       The wavelengths vector in microns
;
;       Lambda0
;       The reference wavelength in microns
;
;

;--------------------------------------------------------
; wavelength range and reference wavelength to edit
;--------------------------------------------------------

Lambda1=0.3		; smallest wavelength [um]
Lambda2=1.0		; largest wavelength [um]
Lambda_Step=0.05        ; wavelength step [um]
Lambda0=0.5		; reference wavelength [um]

;---------------------------------------------------------
; atmospheric parameters to edit
;---------------------------------------------------------

TC=11.5D			;Temperature [C]
RH=14.5D			;Relative Humidity [%]
P=743.0D			;Pressure [mbar]

;----------------------------------------------------------


Lambda=Findgen((Lambda2-Lambda1)/Lambda_Step+1)*Lambda_Step+Lambda1

ZD=Z*!DtoR
T=TC+273.16D

PS=-10474.0+116.43*T-0.43284*T^2+0.00053840*T^3

P2=RH/100.0*PS
P1=P-P2

D1=P1/T*(1.0+P1*(57.90D*1.0E-8-(9.3250D*1.0E-4/T)+(0.25844D/T^2)))
D2=P2/T*(1.0+P2*(1.0+3.7E-4*P2)*(-2.37321E-3+(2.23366/T)-(710.792/T^2)+(7.75141E4/T^3)))

S0=1.0/Lambda0 & S=1.0/Lambda

N0_1=1.0E-8*((2371.34+683939.7/(130-S0^2)+4547.3/(38.9-S0^2))*D1+$
     (6487.31+58.058*S0^2-0.71150*S0^4+0.08851*S0^6)*D2)

N_1=1.0E-8*((2371.34+683939.7/(130-S^2)+4547.3/(38.9-S^2))*D1+$
     (6487.31+58.058*S^2-0.71150*S^4+0.08851*S^6)*D2)

DR=Tan(ZD)*(N0_1-N_1)*206264.8

Return
End

;---------------------------------------------------------------------------

Pro Diff_Atm_Refr

; The routine generates the Differential Atmospheric Refraction table
; using the routine "dar.pro".
; The extrema and the accuracy step of the airmass can be changed
; accordingly to the user purposes.
;
; CALLING SEQUENCE
;
;       Diff_Atm_Refr
;
;
; ARGUMENTS
;
;       None
;
;
; OUTPUT
;
;       The file "diff_atm_refr.dat" containing the refraction data
;       is written on the disk
;

;-------------------------------------------------------------------
; airmass range to edit
;-------------------------------------------------------------------

Air_Mass_Min=1.0	; smallest airmass  >= 1.0 !!! 
Air_Mass_Max=4.9	; largest airmass
Air_Mass_Step=0.05	; airmass step

DAR,0,DR,Lambda,Lambda0
N_Refr=N_Elements(DR)

Lambda=Fix(Lambda*1.0E4)
Lambda0=StrTrim(String(Fix(Lambda0*1.0E4)),1)

WL='      '
For I=0,N_Refr-1 Do Begin
   WL1=StrTrim(String(Lambda(I)),1)
   IF StrLen(WL1) EQ 4 Then WL=WL+'  '+WL1 Else WL=WL+' '+WL1
EndFor

Air_Mass=Findgen((Air_Mass_Max-Air_Mass_Min)/Air_Mass_Step+1)*Air_Mass_Step+Air_Mass_Min
N_AMass=N_Elements(Air_Mass)
Coeffs=FltArr(N_AMass,N_Refr)

Z=ACos(1.0/Air_Mass)*!Radeg

For I=0,N_AMass-1 Do Begin
   DAR,Z(I),DR
   Coeffs(I,*)=-DR
EndFor

OpenW,UnitW,'diff_atm_refr.dat',/Get_LUn

PrintF,UnitW,'Differential atmospheric refraction at Paranal as a function of wavelength'
PrintF,UnitW,'(horizontal, in Angstrom) and airmass (vertical).'
PrintF,UnitW,'The values are in arcseconds with respect to a wavelength of '+Lambda0+' Angstrom' 
PrintF,UnitW,' '
PrintF,UnitW,WL
PrintF,UnitW,' Sec Z'
PrintF,UnitW,' '

   For I=0,N_AMass-1 Do Begin
      PrintF,UnitW,[Air_Mass(I),Reform(Coeffs(I,*),15)],Format='(16(1X,F5.2))'
   EndFor

Free_LUn,UnitW

Return
End