;+
; NAME:
;      inms_ion_transmission - computes the effect of angle on transmission
;
function inms_ion_transmission, nMass,  nVelX, nVelY, nVelZ, debug=debug
;
; ARGUMENTS:
;   nMass:  The atomic mass of the ion species in DA
;   nVelX:  The component of the ion velocity in the boresight direction
;   nVelY:  The component of the velocity in the Y axis direction
;   nVelZ:  The component of the velocity in the Z axis direction
;
; KEYWORDS:
;   /debug  Set to enable debugging
;
; RETURN VALUE:
;   The reduction in transmission due to angle of attack,
;    T(a) = F * T(0)
;   T is transmission, F is the value returned.
;   In the event of an error a value of -1 is returned
;
; CHANGES:
;   28 July 2006  D.Gell Initial coding based on W. Kasprzak routine
;                         Open_Source_Ion_Angle
;    20-Sep-2006  DAG    Use error handler and inms_post_message
;
; CONSTANTS:

    nConvert = 5.1824E-03         ;; kenetic energy conversion

    nZoffset = 0.0                ;; z direction response model
    nZ0 = 3.873E-04
    nZ1 = 4.218E-01

    nYoffset = 0.0                ;; y direction response model
    nY0 = 2.534E-03
    nY1 = 1.049E-01
;-
;===========================================================================
;
    ;; establish error handler
    bDebugSw = keyword_set(debug)
    nError = 0
    catch, nError
    if nError ne 0 then begin
        catch, /cancel
        inms_post_message, traceback=bDebugSw, $
                           console=inms_is_image() or bDebugSW
        if bDebugSw then stop
        return, -1
    endif 

    ;;validate arguments

    if n_params() ne 4 then begin
        message, 'One or more required parameters are missing'
    endif 

    ;; compute speed & velocity projected on y-z plane
    nSpeed = sqrt(nVelX^2 + nVelY^2 +nVelZ^2)

    ;; compute particle energy
    nEnergy = nConvert * nMass * nSpeed^2

    ;; compute angular projections of velocity onto
	;; entrance plane

    nAng1  = !radeg * acos(-nVelX / nSpeed)
    nAng2  = atan(nVelZ, nVelY)

    ;; compute response
    nZwide = nZ0 * exp(nZ1*nEnergy)
    nZresp = exp(-(abs(nAng1-nZoffset)^4.0) * nZwide)

    nYwide = nY0 * exp(nY1*nEnergy)
    nYresp = exp(-(abs(nAng1-nYoffset)^4.0) * nYwide)
  
    anZero = where(nZresp eq 0. or nYresp eq 0., nZero, $
                   complement=anNonZero, ncomplement=nNonZero)
    nResponse = fltarr(n_elements(nSpeed))

    if nZero gt 0 then nResponse[anZero] = 0

    if nNonZero gt 0 then begin
        nResponse[anNonZero] = sqrt((nYresp[anNonZero]*nZresp[anNonZero])^2 $
                                    / ((nZresp[anNonZero]*Cos(nAng2[anNonZero]))^2 $
                                    +  (nYresp[anNonZero]*Sin(nAng2[anNonZero]))^2))
    endif 
    return, nResponse
end