; NAME:
;   add_cosmic_rays
;
; PURPOSE:
;   Simulates cosmic ray hits and returns a 2D array of hits on an array of
;   detector pixels over a certain exposure time.
; 
; INPUTS:
;   xsize       number of detector pixels in x-dimension
;   ysize       number of detector pixels in y-dimension
;   time        time for which detector is exposed to cosmic rays (in seconds)
;   pix         pixdim struct containg x,y,z dimensions of detector pixels
;  
; KEYWORDS:
;   rate        incidence rate for cosmic rays, in #/cm^2/second.  Default value
;               is 5.
;   mask        when set, keyword indicates that charge diffusion to neighboring
;               pixels should be modeled using a hard-coded diffusion pattern
;               mask.
;
; OUTPUTS:
;   crimage     floating point array of cosmic ray hits
;
; REFERENCE:
;   This code is adapted from a method in the multi-data sim package by 
;   J. Offenberg
;
Function add_cosmic_rays, xsize, ysize, time, pix, rate = rate, mask = mask

    ; A static random number seed is needed, but IDL doesn't have static
    ; variables.  A common block can be used to get the same effect--the value
    ; of the crseed variable will persist over different executions of this
    ; function.
    common crstatic, crseed

    ; if the rate is not specified, use the default of 5 rays/cm^2/sec.
    if (not(keyword_set(rate))) then begin
        rate = 5.0
    endif
        
    ; calculate mean number of rays hitting detector in time interval
    n_rays = xsize * (pix.x * 1e-4) * ysize * (pix.y * 1e-4) * rate * time

    ; Replace the mean n_rays with a value from a Poisson distribution about 
    ; the mean to account for shot noise.
    n_rays = randomu(crseed, poisson = n_rays)

    ; Create an empty cosmic ray image
    crimage = fltarr(xsize, ysize)

    ; define the differential unit of path length of the ray through the
    ; detector (in microns)
    dpath = 0.1

    ; define cosmic ray mask used to model charge diffusion
    if (keyword_set(mask)) then begin
        crmask = [ $
            [1.5e-3, 1.6e-2, 1.5e-3], $
            [1.6e-2, 1.0000, 1.6e-2], $
            [1.5e-3, 1.6e-2, 1.5e-3]]
        poimask = fltarr(3,3)
    endif
    
    ; Generate the rays
    for i = 0, n_rays - 1 do begin

        ; calculate random impact point for the ray
        xpos = randomu(crseed) * xsize
        ypos = randomu(crseed) * ysize
        zpos = 0.0

        ; calculate angles of incidence
        theta = randomu(crseed) * 2.0 * !pi
        cosphi = randomu(crseed)
        sinphi = sqrt(1.0 - cosphi^2)
        phi = asin(sinphi)

        ; define the mean number of electrons liberated for each dpath length
        ; that the ray travels through the detector.
        ; 10% of the time the cosmic ray is a He nucleus, with 4 times the
        ; effect of a proton.
        if (randomu(crseed) ge 0.9) then begin
            depath = 400.0
        endif else begin
            depath = 100.0
        endelse

        ; calculate x, y and z components of dpath
        dx = cos(theta) * dpath * sinphi
        dy = sin(theta) * dpath * sinphi
        dz = dpath * cosphi


        ; compute position of cosmic ray at each step and add generated signal
        ; to detector
        while (zpos le pix.z and xpos ge 0.0 and ypos ge 0.0 and xpos lt xsize $
            and ypos lt ysize) do begin
  

            if (keyword_set(mask)) then begin

                x0 = (fix(xpos) - 1) > 0
                x1 = (fix(xpos) + 1) < (xsize - 1)
                y0 = (fix(ypos) - 1) > 0
                y1 = (fix(ypos) + 1) < (ysize - 1)

                mx0 = x0 - (fix(xpos) - 1)
                mx1 = mx0 + (x1 - x0)
                my0 = y0 - (fix(ypos) - 1)
                my1 = my0 + (y1 - y0)
                

                for mx = mx0, mx1 do begin
                    for my = my0, my1 do begin
                        poimask[mx, my] = randomu(crseed, $
                            poisson = depath * crmask[mx, my])
                    endfor
                endfor

                crimage[x0:x1, y0:y1] = crimage[x0:x1, y0:y1] + $
                    poimask[mx0:mx1, my0:my1]

            endif else begin

                crimage[fix(xpos), fix(ypos)] = crimage[fix(xpos),fix(ypos)] + $
                    randomu(crseed, poisson = depath)

            endelse

            xpos = xpos + (dx / pix.x)
            ypos = ypos + (dy / pix.y)
            zpos = zpos + dz

        endwhile        

        ; account for the last sub-interval that must be traversed to reach the
        ; bottom of the detector surface.
        zrem = (pix.z - zpos + dz)
        if ((xpos ge 0) and (ypos ge 0) and (xpos lt xsize) and $
            (ypos lt ysize) and (zrem gt 0.0)) then begin
            depath = depath * zrem / dz
            crimage[fix(xpos), fix(ypos)] = crimage[fix(xpos), fix(ypos)] + $
                randomu(crseed, poisson = depath)
        endif


    endfor


    return, crimage
        
End