pro set_result,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 CHAP,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.E20
  A1 = 1.0606963
  A2 = .55643831
  A3 = 1.0619896
  A4 = 1.7245609
  B1 = .56498823
  B2 = .06651874

;  Set up result
  result=(theta*0.)-999.
;  vectorize scale height and altitude if necessary
  h=h_+(theta*0)
  z=z_+(theta*0)
;  Check the calling arguments for reasonableness
  bad=logical_or(logical_or(Z LT 0.,H LT 0.),ABS(THETA) GT 180.)
  set_result,result,bad,0
;
;  Calculate the value
  set_result,result,theta eq 0,1
  set_result,result,THETA GT (90. + SQRT(Z)),CHAP_MIN
  thetabig=abs(theta) gt 90
  w_thetabig=where(thetabig,n_thetabig,complement=w_nthetabig,ncomplement=n_nthetabig)
  angle=result*0.;
  F=result*0.;
  temp=result*0.;
  IF n_thetabig gt 0 THEN BEGIN
    ANGLE[w_thetabig] = !PI - ABS(THETA[w_thetabig]*!const.dtor)
    F[w_thetabig] = -1.
    rp=result*0.
    ARG=result*0.
    RP[w_thetabig] = (RE+Z[w_thetabig])*SIN(ANGLE[w_thetabig])
    ARG[w_thetabig] = (RE+Z[w_thetabig]-RP[w_thetabig])/H[w_thetabig]
    set_result,result,arg gt 170,chap_min
    TEMP[w_thetabig] = SQRT(2.*!PI*RP[w_thetabig]/H[w_thetabig])*EXP(ARG[w_thetabig])
  end
  if n_nthetabig gt 0 then BEGIN
    F[w_nthetabig] = 1.
    ANGLE[w_nthetabig] = THETA[w_nthetabig]*!const.dtor
    TEMP[w_nthetabig] = 0.
  END
  C = SQRT(!PI*(RE+Z)/(2.*H))
  Y = C*COS(ANGLE)/SQRT(!PI)
  set_result,result,y lt 8,TEMP + F*C*(A1 + A2*Y)/(A3 + (A4 + Y)*Y)
  set_result,result,~(y lt 8),TEMP + F*C*B1/(B2 + Y)
  return,result
END