Commit 752befc7 authored by Thierry Corbard's avatar Thierry Corbard

import public files from svn

parent f6c85ff9
;+
; NAME: air_index
;
; PURPOSE:
; Compute ambiant air refractive index following
; either Baldini (1963) (IUGG 1963 adopted standard) or Ciddor (1996) (IAG
; 1999 adopted standard)
;
; CATEGORY:
; Atmospheric Physics
;
; CALLING SEQUENCE:
; nobs=air_index(Ta,P,hum,lbda,type[,xc=xc])
;
; INPUTS:
; Ta: [K] temperature
; P: [Pa] Pressure
; hum: humidity [%] in [0,100]
; lbda: wavelength [m]
; type: 'Baldini' or 'IUGG1963' <br>
; 'Ciddor' or 'IAG1999'
;
; OPTIONAL INPUTS:
;
; KEYWORD PARAMETERS:
; xc: [ppm] Carbon dioxyde concentration (used in 'Ciddor' model only). <br>
; 380 ppm is taken by default (if xc is not set) in 'Ciddor' model. <br>
; 300 ppm is assumed in 'Baldani' model.
;
; OUTPUTS:
; Refractive index $n$ of ambiant air at given wavelength.
;
; PROCEDURE:
;
; +Model 'Baldini'/'IUGG1963' (IUGG, 1963, depreciated) <br><br>
;
;\begin{equation}
;T_0=273.15\, \mathrm{K}\, ,P_0=101325\, \mathrm{Pa}\, ,x_c=300\, \mathrm{ppm}
;\end{equation}
; \begin{equation}
; n(T,P,f_h,\lambda)-1={T_0 \over T} \left\{(n_0(\lambda)-1) {P
; \over P_0} - 4.13\, 10^{-10}\ p(f_h,T)\right\} \,\, [1, 2]
; \end{equation}
;\begin{equation}
;n_0(\lambda)-1=\left\{2876.04+{16.288 \over
;(10^6\lambda)^2}+{0.136\over (10^6\lambda)^4}\right\}\ 10^{-7} \,\,[3]
;\end{equation}
;\begin{equation}
;p(f_h,T)=f_h\ 6.1075\, 10^2\ e^{7.292\,
;10^{-2}(T-T_0)-2.84\,10^{-4}(T-T_0)^2} \,\, [4,5]
;\end{equation}
; [1] Baldini, A. A. 1963, GIMRADA Research Note, 8 <br>
; [2] IUGG (International Union of Geodesy and Geophysics), 1963, Resolutions,
; 13th General Assembly, 19-31 August 1963, Berkeley, California, USA.
; Bulletin Geodesique, 70, 390 <br>
; [3] Barrel, H. & Sears, J. E. 1939, Phil. Trans. R. Soc., A238, p. 52,
; equation (7.7) for f=0, t=0, p=760 mmHg
; (assumed xc=300 ppm (0.03%) ) <br>
; [4] Chollet F. 1981, PhD thesis, Universite Pierre et Marie Curie (Paris VI) <br>
; [5] Annuaire du Bureau Des Longitudes, 1975 <br><br>
; +Model 'Ciddor'/'IAG1999' (IAG adopted standard, 1999) <br><br>
; [6] Ciddor, P.E., 1996, "Refractive index of air: new equations for the visible and near infrared", Applied Optics LP, vol. 35, Issue 9, p.1566 <br>
; [7] <a href="http://iag.dgfi.tum.de/fileadmin/IAG-docs/IAG_Resolutions_1999.pdf" > IAG
; (International Association of Geodesy), 1999, Resolution 3, 22nd
; GeneralAssembly, 19-30 July 1999,Birmingham, U.K.</a>
; <br> CALLS:
; $SSW/optical/gen/idl/util/refractivity.pro (SolarSoft library)
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
; Written by: T. Corbard 06/2012 (corbard@oca.eu)
;
;-
function air_index,Ta,P,hum,lbda,type,xc=xc
COMPILE_OPT IDL2
IF (N_ELEMENTS(xc) NE 1) THEN xc = 380.d0
T0=273.15D0
t=Ta-T0
if(type eq 'Cidor') then type='Ciddor'
case 1 of
(type eq 'Baldini') or (type eq 'IUGG1963'): BEGIN
P0=101325.D0 ;Pa
lbd=lbda*1.d6 ;micrometer
;al=1.D0/273.15d0;0.00367d0
lbd2=lbd*lbd
lbd4=lbd2*lbd2
n0m1=1.D-7*(2876.04D0+16.288D0/lbd2+0.136D0/lbd4) ;[3]
hsat=4.581D0*exp(7.292D-2*t-2.84D-4*t^2) ;mmHg [4, 5]
eh=(hum*hsat/100.D0)
;n=1.D0+((1.D0/(al*Ta))*(n0m1*(P/P0)-5.5D-8*eh))
n=1.D0+((T0/Ta)*(n0m1*(P/P0)-5.5D-8*eh)) ;[1, 2]
END
(type eq 'Ciddor') or (type eq 'IAG1999'):BEGIN
wavelength=lbda*1.d9
nm1=refractivity(wavelength, t, P, hum,xc) ;[6,7]
n=nm1+1.D0
END
ELSE: message,'Unknown type: '+type
ENDCASE
return,n
end
function curvature_radius,lat,h,Az=Az
;+
; NAME: curvature_radius
;
; PURPOSE: Compute various estimate of the radius of curvature
; and distance from geocenter at observer position
;
; CATEGORY:
; Geophysics
;
; CALLING SEQUENCE:
; radius=curvature_radius(lat,h[,Az=Az])
;
; INPUTS:
; lat: Geodetic latitude at observer [rad]
; h: Elevation of observer above the reference ellipsoid (added to the output) [m]
;
; OPTIONAL INPUTS:
;
; KEYWORD PARAMETERS:
; Az: Azimuth angle [rad] (default =0)
;
; OUTPUTS:
; radius[0] curvature in the (noth-south) meridian [m] <br>
; radius[1] curvature in the prime vertical [m] <br>
; radius[3] mean radius of curvature (=sqrt(radius[0]*radius[1])) [m] <br>
; radius[4] distance from geocenter (WGS84) [m] <br>
; radius[5] radius of curvature at azimuth angle Az (if set) [m]
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
; Use WGS84 Reference ellipsoid a=6378137 m, b=6356752 m <br>
;\begin{equation}r_c^0={(ab)^2 \over
;{\left(a^2\cos^2(\varphi)+b^2\sin^2(\varphi) \right)^{3/2}}}
;\end{equation}<br><br>
;
;\begin{equation}
;r_c^{90}={a^2 \over {\sqrt{a^2\cos^2(\varphi)+b^2\sin^2(\varphi)} }}
;\end{equation}<br><br>
;
; \begin{equation}
;<\!\!r_c\!\!>=\sqrt{r_c^0 r_c^{90}}= {a^2b \over
;{a^2\cos^2(\varphi)+b^2\sin^2(\varphi) }}
;\end{equation}<br><br>
;
;\begin{equation}
;R=\sqrt{ {{a^4\cos^2(\varphi)+b^4\sin^2(\varphi)} }\over
;{{a^2\cos^2(\varphi)+b^2\sin^2(\varphi)} }}
;\end{equation} <br><br>
;
;\begin{equation}
;r_c^A={1 \over {{\cos^2(A)\over r_c^0}+{\sin^2(A)\over r_c^{90}} } }
;\end{equation}
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
; T. Corbard 06/2012 corbard@oca.eu
;
;-
COMPILE_OPT IDL2
;WGS84 reference ellipsoid
a=6378137.D0 ;[m]
b=6356752.D0 ;[m]
a2=a*a
b2=b*b
sl2=sin(lat)^2
cl2=cos(lat)^2
c=a2*cl2+b2*sl2
r0=a2*b2/c^(1.5)
r90=a2/sqrt(c)
rm=sqrt(r0*r90)
Rd=sqrt((a2*a2*cl2+b2*b2*sl2)/c)
if (N_elements(Az) ne 0) then begin
rA=1.D0/(cos(Az)^2/r0+sin(Az)^2/r90)
return,[r0,r90,rm,Rd,rA]+h
endif else begin
return,[r0,r90,rm,Rd]+h
endelse
end
function meteo,model=model,T=T,P=P,hum=hum,laps=laps,xc=xc,density=density
;+
; NAME:
; meteo
;
; PURPOSE:
; gives a structure that defines atmospheric conditions
;
; CATEGORY:
; Atmospheric Physics
;
;
; KEYWORD PARAMETERS:
; model: string for predefined weather conditions <br>
; 'Std' T=15$^\circ$C, P=1013.25 hPa, hum=0 %, laps=6.5 K/Km, xc=450 ppm <br>
; 'Isothermal' T=15$^\circ$C, P=1013.25 hPa, hum=0 %, laps=1.d-12 K/m, xc=450 ppm <br>
; 'Adiabatic' T=15$^\circ$C, P=1013.25 hPa, hum=0 %, laps=10 K/km, xc=450 ppm <br>
; 'Calern_mean' T=15$^\circ$C, P=875 hPa, hum=50 %, laps=6.5 K/km, xc=380 ppm <br>
; 'Calern_min' T=-10$^\circ$C, P=900 hPa, hum=50 %, laps=6.5 K/km, xc=380 ppm <br>
; 'Calern_max' T=30$^\circ$C, P=850 hPa, hum=50 %, laps=6.5 K/km, xc=380 ppm <br><br>
; T: [K] Temperature
; P: [Pa] Pressure
; hum: [%] in [0-100] humidity
; laps: [K/m] Troposheric laps rate
; xc: [ppm] Carbon dioxid concentration
; density: [kg/m^3] Air density
;
; OUTPUTS:
; Structure containing temperature (meteo.T), pressure (meteo.P),
; humidity (meteo.hum), troposheric
; laps rate (meteo.laps) and CO2 concentration (meteo.xc). <br>
; Missing information from input are replaced by standard atmosphere
; values (Ciddor 1996) or by the chosen model if [model=] is set.
;
; EXAMPLE:
; std_atm=meteo() => return Ciddor 1996 standard atmosphere parameters.
; atm1=meteo(model='Isothermal',T=293.D0,P=100000.D0) Isothermal model
; except for T and P.
;
; MODIFICATION HISTORY:
; T. Corbard 02/2013 corbard@oca.eu
;
;-
;Standard Atmosphere in Ciddor (1996)
std_atm ={meteo,T:273.15D0+15.D0,P:101325.D0,hum:0.D0,laps:0.0065D0,xc:450.D0, density:0.D0}
isothermal ={meteo,T:273.15D0+15.D0,P:101325.D0,hum:0.D0,laps:1.D-12 ,xc:450.D0, density:0.D0}
adiabatic= {meteo,T:273.15D0+15.D0,P:101325.D0,hum:0.D0,laps:0.010D0,xc:450.D0, density:0.D0}
Calern_mean={meteo,T:273.15D0+15.D0,P:87500.D0,hum:50.D0,laps:0.0065D0,xc:380.D0, density:0.D0}
Calern_min= {meteo,T:273.15D0-10.D0,P:90000.D0,hum:50.D0,laps:0.0065D0,xc:380.D0, density:0.D0}
Calern_max= {meteo,T:273.15D0+30.D0,P:85000.D0,hum:50.D0,laps:0.0065D0,xc:380.D0, density:0.D0}
if(N_elements(model) eq 0) then model='Std'
CASE model OF
'Std':atm=std_atm
'Isothermal':atm=isothermal
'Adiabatic':atm=adiabatic
'Calern_mean':atm=Calern_mean
'Calern_min':atm=Calern_min
'Calern_max':atm=Calern_max
ELSE : message, 'unknown meteo: '+model
ENDCASE
if(N_elements(T) ne 0) then if(finite(T)) then atm.T=T
if(N_elements(P) ne 0)then if(finite(P)) then atm.P=P
if(N_elements(hum) ne 0) then if(finite(hum)) then atm.hum=hum
if(N_elements(xc) ne 0) then if(finite(xc)) then atm.xc=xc
if(N_elements(laps) ne 0)then if(finite(laps)) then atm.laps=laps
atm.density=atmospheric_density(atm.T-273.15D0, atm.P, atm.hum, atm.xc)
if(N_elements(density) ne 0)then if(finite(density)) then atm.density=density
return,atm
end
function reduced_height,arg
;+
; NAME:
; reduced_height
;
; PURPOSE:
; Give the "reduced height" or height of an homogenous atmosphere
; with the same density as the density at observer position and that
; would give the same pressure as the one recorded at observer position.
;
; CATEGORY:
; Atmospheric Physics
;
; CALLING SEQUENCE:
; l=reduced_height([T])
; l=reduced_height([P,rho,g])
;
; INPUTS:
; arg: 1 or 3 element vector <br>
; arg[0] : Temperature [K] or
; Pressure [Pa] (if arg[1:2] are defined) <br>
; arg[1]: Density [kg/m$^3$] (optional)<br>
; arg[2]: gravity [m/s$^2$] (optional) <br>
;
; OPTIONAL INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
; homogeneous atmosphere reduced height [m]
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
; \begin{equation}
; \ell(P,\rho,g)={P \over {\rho \ g}}
; \end{equation}
; If only temperature is given, ideal gas law for dry air
; is assumed.
; \begin{equation}
;\ell(T)={P_0 \over {\rho_0 \ g_0}}{T \over T_0}
;\end{equation}
;where ${\rho_0=1.293\ \mbox{kg}\,\mbox{m}^{-3}}$ for
;${T_0=273.15\,\mbox{K}}$, ${P_0=101325\,\mbox{Pa}}$ and normal
;gravity ${g_0=9.80665\,\mbox{m}\,\mbox{s}^{-2}}$.
;
; EXAMPLE:
; ;Compare reduced height assuming ideal gaz law or from full
; ;air density calculation
;
;t=15.D0 ;Celsius
;P=875d2 ;Pa
;hum=50.D0 ;%
;xc=380 ; ppm
;g=9.802155492d0 ; m/s^2 Measured gravity acceleration at Calern
; ;Get the air density from Ciddor (1996) equations
;d=atmospheric_density(t,P,hum,xc)
;print,'Air Density: ',d,' kg/m^2',format='(a,d5.3,a)'
;print,'Reduced height (Full density calculation): ',round(reduced_height([P,d,g])),' m'
;print,'Reduced height (Assuming ideal gaz law ): ',round(reduced_height([t+273.15D0])),' m'
;
;Air Density: 1.054 kg/m^2
;Reduced height (full density calculation): 8467 m
;Reduced height (assuming ideal gaz law ): 8430 m
;
; MODIFICATION HISTORY:
; T. Corbard 06/2012 corbard@oca.eu
;
;-
l=0
CASE N_elements(arg) of
1: BEGIN
Ta=arg[0] ;K
T0=273.15d0 ;K
P0=101325.D0 ;Pa
g0=9.80665 ;m s^-2
rho0=1.293 ;kg m^-3
l=P0/(rho0*g0)*Ta/T0 ;m
END
3: BEGIN
P=arg[0] ;Pa
rho=arg[1] ;kg m^-3
g=arg[2] ;m s^-1
l=P/(rho*g) ;m
END
else: message,'unvalid argument'
ENDCASE
return,l
end
This diff is collapsed.
function refrac_model,name
;+
; NAME:
; refrac_model
;
; PURPOSE:
; returns a structure that defines a refraction model
;
; CATEGORY:
; Astronomy
;
; CALLING SEQUENCE:
; refraction=refrac_model('Model_name')
;
; INPUTS:
; name: string giving the refraction model name (case insensitive).
; Available models are : <br>
; 'NAO63Ciddor' <br>
; 'tan5Baldini' <br>
; 'tanCiddor' <br>
; 'tan5Ciddor' <br>
; 'erfcCiddor' <br>
; 'CassiniCiddor' <br>
; 'Laplace' (same as 'tan3Ciddor') <br>
; 'Laclare' (same as 'tan5Baldini') <br><br><br>
;
; the first chararters give the type of refraction model: <br><br>
;
; 'NAO63' <br>Exact evaluation of the refraction integral from
; standard atmospheric model.
; N.A.O Technical Notes 59 and 63 and a paper by Standish and Auer
; "Astronomical Refraction: Computational Method for
; all Zenith Angles". Modified to take Cidor (1996) modern
; air index calculation (see <a href="refrac_all_p.html#aref" >aref.pro </a> ) <br><br>
;
; 'Laclare' <br>Code provided by F. Laclare for Calern Observatory
; (see <a href="refrac_all_p.html#refrac_laclare" >refrac_laclare.pro </a> )
; <br><br>
; or one of the formulae involving only the atmospheric air index and
; reduced height (see <a href="refrac_all_p.html#refrac_for" >refrac_for.pro </a> ) :<br><br>
; 'tan' One term expansion R=k*tan(z) <br>
;
; 'tan3' Laplace Formula (2 terms expansion in odd powers of
; tan(z)) <br>
;
; 'tan5' 3 terms expansion (Danjon 1980, derived from erfc formula) <br>
;
; 'erfc' (error function, Danjon 1980, assuming exponential
; law for variation of density with height) <br>
;
; 'Cassini' Cassini's model for homogeneous atmosphere <br><br>
;
; The following characters give the type of model for air index (see <a href="air_index.html#air_index" >air_index.pro </a> ) (not used if type='NAO63' or
; type='Laclare') <br><br>
;
; 'Baldani' (Chollet 1991 with Barel and Sears (1939) for
; refractivity and IUGG 1963 +BDL 75 for pressure,
; temperature, humidity dependence) <br><br>
;
; 'Ciddor' Ciddor (1996) modern calculation <br><br>
;
; OUTPUTS:
; structure defining the refraction model
; (input for procedure refrac_all_p.pro)
;
; EXAMPLE:
; ref_NAO63Ciddor=refrac_model('NAO63Ciddor')
; ref_Laclare=refrac_model('Laclare')
;
; MODIFICATION HISTORY:
; T. Corbard 02/2013 corbard@oca.eu
;
;-
;Check if 'Cidor' is used in place of 'Ciddor'
namelow=STRLOWCASE(name)
p=strpos(namelow,'cidor')
if(p ne -1) then namelow=strmid(namelow,0,p)+'ciddor'+strmid(namelow,p+5,strlen(namelow)-p-4)
CASE namelow OF
'tan3ciddor' :ref={refrac, type:'tan3' ,type_rc:'meridian',type_ai:'Ciddor' ,Az:0.D0}
'laplace' :ref={refrac, type:'tan3' ,type_rc:'meridian',type_ai:'Ciddor' ,Az:0.D0}
'laclare' :ref={refrac, type:'Laclare',type_rc:'' ,type_ai:'' ,Az:0.D0}
'nao63ciddor' :ref={refrac, type:'NAO63' ,type_rc:'meridian',type_ai:'' ,Az:0.D0}
'tan5baldini' :ref={refrac, type:'tan5' ,type_rc:'meridian',type_ai:'Baldini' ,Az:0.D0}
'tanciddor' :ref={refrac, type:'tan' ,type_rc:'meridian',type_ai:'Ciddor' ,Az:0.D0}
'tan5ciddor' :ref={refrac, type:'tan5' ,type_rc:'meridian',type_ai:'Ciddor' ,Az:0.D0}
'erfcciddor' :ref={refrac, type:'erfc' ,type_rc:'meridian',type_ai:'Ciddor' ,Az:0.D0}
'cassiniciddor':ref={refrac, type:'Cassini',type_rc:'meridian',type_ai:'Ciddor' ,Az:0.D0}
ELSE : message, 'unknown refraction model: '+name
ENDCASE
return,ref
end
function station,name,h=h,lat=lat
;+
; NAME:
; station
;
; PURPOSE:
; define a structure that gives elevation above the reference
; ellipsoid (in meters) and geodetic latitude (in
; radian) of a given observing site
;
; CATEGORY:
; Astronomy
;
; CALLING SEQUENCE:
; site=station('Site_name') ;look in the database
; site=station('Site_name',h=h,lat=[deg,min,sec]) ;just create the
; structure from user inputs h and lat
;
; INPUTS:
; name: string giving the Site name. <br>
;
; KEYWORD PARAMETERS:
; h: altitude above WGS84 reference ellipsoid in meters
; lat: geodetic latitude either in fractional degrees or as a vector [degree,min,sec]
;
; OUTPUTS:
; 'station' Structure <br>
; station.h (elevation above the reference ellipsoid) <br>
; station.lat (geodetic latitude in rad) <br>
; station.name (string) name of the station
;
;PROCEDURE:
;
; define station structure for station in the database or use keywords to create user defined
; station structure<br><br>
;
; CALLS: <a href="http://idlastro.gsfc.nasa.gov/ftp/pro/astro/ten.pro"> ten.pro </a> from IdlAstro lib
;
; EXAMPLE:
;
; calern=station('Calern') ;provide geodetic latitude and elvation above the ellipsoid
; calern2=station('Calern2',h=1271.d0,lat=[43.D0,44.D0,53.D0]) ;astronomic latitude and elevation above see level
; ;Print various radius and show difference in meters
; print,round(curvature_radius(calern.lat, calern.h))
; print,round(curvature_radius(calern2.lat,calern2.h)-curvature_radius(calern.lat, calern.h))
;
; 6367308 6389694 6378491 6369278
; -56 -53 -55 -51
;
; MODIFICATION HISTORY:
; T. Corbard 02/2013 corbard@oca.eu <br>
; TC Oct. 2015 Added call to IDLAstro 'ten'. Allows
; user defined station with optional keywords h and lat <br>
; Replace the astronomical latitude by the geodetic latitude of Calern <br>
; Replace the altitude above sea level by the altitude above the WGS84
; reference ellipsoid for Calern.
;-
if((N_elements(h) ne 0) && (N_elements(lat) ne 0)) then begin
sta={station, h:h, lat:ten(lat)*!DPI/180.D0, name:name}
endif else begin
if((N_elements(h) ne 0) || (N_elements(lat) ne 0)) then message,'Both h and lat must be provided'
CASE STRLOWCASE(name) OF
;'calern':sta={station, h:1274.2423d0,lat:(43.D0+(44.D0+((55.D0)/60.D0))/60.D0)*!DPI/180.D0, name:'Calern'}
'calern':sta={station, h:1323.d0, lat:ten([43.D0,45.D0,7.D0])*!DPI/180.D0, name:'Calern'}
ELSE : message, 'Unknown station, please use keywords h= and lat='
ENDCASE
END
return,sta
END
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment