FUNCTION CHAP_ORIG,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
;
;  Input:
;    Z     - Altitude (km)
;    THETA - Solar zenith angle (degrees)
;    H     - Scale height (km)
;
;  Output:
;    Chapman's grazing incidence integral returned as function value
;
    COMMON CHAP_COM,RE,CHAP_MIN,A1,A2,A3,A4,B1,B2
;
;  Initialize COMMON CHAP_COM if this is the first call
    ASIZE = SIZE(RE)
    IF ASIZE(1) EQ 0 THEN BEGIN
       RE = 6371.
       CHAP_MIN = 1.E-20
       A1 = 1.0606963
       A2 = .55643831
       A3 = 1.0619896
       A4 = 1.7245609
       B1 = .56498823
       B2 = .06651874
    ENDIF
;
;  Check the calling arguments for reasonableness
    IF (Z LT 0.) OR (H LT 0.) OR (ABS(THETA) GT 180.) THEN BEGIN
       PRINT,'Error in arguments to Chapman function'
       PRINT,Z,THETA,H,FORMAT= $
             "('Z = ',F5.1,5X,'Theta = ',F5.1,5X,'Scale height = ',F9.1)"
       RETURN,0.
    ENDIF
;
;  Calculate the value
    IF THETA EQ 0. THEN RETURN,1. $
       ELSE IF THETA GT (90. + SQRT(Z)) THEN RETURN,CHAP_MIN $
       ELSE BEGIN
          IF ABS(THETA) GT 90. THEN BEGIN
             ANGLE = !PI - ABS(THETA*!const.dtor)
             F = -1.
             RP = (RE+Z)*SIN(ANGLE)
             ARG = (RE+Z-RP)/H
             IF ARG GT 170. THEN RETURN,CHAP_MIN
             TEMP = SQRT(2.*!PI*RP/H)*EXP(ARG)
          ENDIF ELSE BEGIN
             F = 1.
             ANGLE = THETA*!const.dtor
             TEMP = 0.
          ENDELSE
          C = SQRT(!PI*(RE+Z)/(2.*H))
          Y = C*COS(ANGLE)/SQRT(!PI)
          IF Y LT 8. THEN CHAP = TEMP + F*C*(A1 + A2*Y)/(A3 + (A4 + Y)*Y) $
                     ELSE CHAP = TEMP + F*C*B1/(B2 + Y)
          RETURN,CHAP
       ENDELSE
;
    END