FUNCTION vanrhijn,zenang,emis_ht,re=re
;+
; NAME:
;	VANRHIJN
; PURPOSE:
;       Generation of a van Rhijn profile for a uniformly bright
;       emitting layer seen at zenith angles given by the zenang
;	array.
; CATEGORY:
;	Image correction.
; CALLING SEQUENCE:
;	vR = VANRHIJN( ZENANG, EMIS_HT, RE=RE)
; INPUTS:
;       ZENANG: A variable (scalar or array) giving the zenith
;               angles at which the sky is viewed.  These must be
;		in DEGREES.
;       EMIS_HT: Scalar variable giving the height of the emitting
;		layer, in km.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;       RE:     The local radius of the Earth.  The default value
;		is 6371.2 km.
; OUTPUTS:
;       vR:     A variable of the same structure as ZENANG, giving
;               the column-integrated emission brightnesses, at
;               zenith angles ZENANG, of a uniformly-emitting layer
;		of unit brightness when seen in the zenith.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;       If zenith angles greater than 89 degrees are passed, they
;	are truncated to 89 deg.
; PROCEDURE:
;	Straightforward.
; EXAMPLE:
;	; Generate an elevation angle profile from 5 to 175 deg elevation.
;	elev=FINDGEN(171)+5.
;	; Convert to zenith angles.
;	za=ABS(90.-elev)
;       ; Compute van Rhijn profile over the meridian for an emitting
;	; layer at 97 km.
;	vr=vanRhijn(za,97.)
; SEE ALSO:
;
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, May 1994.
;	Documented:	D P Steele, July 1994.
;-

; This little function takes zenith angle in degrees (scalar or
; array) and emission height (km) and returns the factor by which
; the van Rhijn effect would increase the apparent column emission
; rate at the given zenith angle(s) over the vertical column emission
; rates.  The calculation is based on the assumption of a thin,
; horizontally uniform emitting layer centered at the emission height.
; Unless the user specifies an Earth radius with the 're' keyword,
; a mean radius of 6371.2 km is assumed.

a=6371.2D0
IF KEYWORD_SET(re) THEN a=DOUBLE(re)

esize=SIZE(emis_ht)
IF esize(0) GT 0 THEN BEGIN
    MESSAGE,'Emission height argument must be a scalar',/INFORMATIONAL
    RETURN,-1
    ENDIF

zsize=SIZE(zenang)
scalar=(zsize(0) EQ 0)
IF scalar THEN zenang=[zenang]

za=zenang<89.

t=(a/(a+emis_ht))^2
v=1./SQRT(1.0D0-t*SIN(!DTOR*za)^2)
IF scalar THEN v=v(0)

RETURN,FLOAT(v)
END