pro set_result2,result,cond,set
  w_set=where(logical_and(result eq -999.,cond),n_set)
  if n_set gt 0 then result[w_set]=set[w_set]
end

FUNCTION CHAP2,Z,THETA,H
;
; Routine to calculate the Chapman function
;  Author: Dave Rusch
;  Modified by Lou Clark, June 1982
;  Converted to IDL by Michael T. Callan, September 2005
;  Vectorized by Chris Jeppesen, February 2006
;
;  Input:
;    Z     - Altitude (km)
;    THETA - Solar zenith angle (degrees)
;    H     - Scale height (km)
;
;  Output:
;    Chapman's grazing incidence integral returned as function value
;
;
;  Initialize COMMON CHAP_COM if this is the first call
  RE = 6371.
  CHAP_MIN = 1.E-20
  A1 = 1.0606963
  A2 = .55643831
  A3 = 1.0619896
  A4 = 1.7245609
  B1 = .56498823
  B2 = .06651874

;  Set up result
  result=(theta*0.)-999.
;
;  Check the calling arguments for reasonableness

  bad=logical_or(logical_or(logical_or(Z LT 0.,H LT 0.),THETA GT 180.),theta lt 0)
  set_result,result,bad,!values.f_nan
;
;  Calculate the value
  set_result,result,theta eq 0,1
  bigmask=logical_and(~bad, theta gt 90)
  w_big=where(bigmask,n_big)
  theta_use=theta
  if n_big gt 0 then theta_use[w_big]=180-theta_use[w_big]
  F = 1.
  ANGLE = THETA_use*!const.dtor
  TEMP = 0.
  C = SQRT(!PI*(RE+Z)/(2.*H))
  Y = C*COS(ANGLE)/SQRT(!PI)
  set_result2,result,y lt 8,TEMP + F*C*(A1 + A2*Y)/(A3 + (A4 + Y)*Y)
  set_result2,result,y ge 8,TEMP + F*C*B1/(B2 + Y)
  if n_big gt 0 then result[w_big]=2*chap(z,90,h)-result[w_big]
  RETURN,result
END