;+
; NAME:
;  fresneldiff_edge_func
; PURPOSE: (one line)
;  Calculate Fresnel diffraction by an edge for curvefit
; DESCRIPTION:
;  Calculate Fresnel diffraction by an edge (opaque for x<0, clear for x>0)
; CATEGORY:
;  Math
; CALLING SEQUENCE:
;  fresneldiff_edge,in,a,f,pder
; INPUTS:
;  in - input structure
;    in.tmid - midtimes (eg seconds)
;    in.dt   - bin widths (seconds)
;    in.tf   - Fresnel width, seconds = sqrt(lam*dist/2)/v_perp
;             since v_perp = dr/dt, this is
;             positive for emersion, negative for immersion
;  a - parameters
;    a[0]   - edge time
;   OR   
;    a[0]   - edge time
;    a[1]   - scale = upper baseline - lower baseline 
;    a[2]   - offset = lower baseline
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  f - Fresnel diffraction pattern, integrated over the bins
;  pder - derivative matrix
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
;  See Goodman Introduction to Fourier Optics, chapter 4
; MODIFICATION HISTORY:
;  Written 2006 Jan 10
;-

pro fresneldiff_edge_func,in,a,f,pder

    ; pull out the description of the problem
    tmid = in.tmid  ; midtimes of bins
    dt = in.dt      ; width of bins
    tf = in.tf      ; fresnel width 
    n = n_elements(tmid)

    deriv = (n_params() gt 3)

    te = a[0] ; time of the edge
    inten = dblarr(n)  ; array for intensity
    if deriv then dinten = dblarr(n) ; array for d(intensity)/d(te)

;    print, 'te = ', te
    ; loop over the bins
    for i = 0L, n-1 do begin
        ; get the edges of this bin
        t0 = tmid[i] - dt[i]/2.d ; start, in sec 
        t1 = tmid[i] + dt[i]/2.d ; end, in sec
        x0 = (t0-te)/tf  ; start, in Fresnel lengths
        x1 = (t1-te)/tf  ; end, in Fresnel lengths

        ; The Fresnel integrals C and S have minima and maxima at 
        ; integral values of x^2
        n_cycles = abs(sgn(x1)*x1^2.d - sgn(x0)*x0^2.)
        if n_cycles gt 100 and sgn(x1) eq sgn(x0) then begin
            inten[i] = 0.5 + 0.5*sgn(x1)
            if deriv then dinten[i] = 0.
        endif else begin
            m = ceil(n_cycles * 3.) ; 3 points per quarter cycle
            dx = (x1-x0)/m
            x = dindgen(m) * dx + x0 + dx/2.d
            c = Fresnel_int(double(x))
            s = Fresnel_int(double(x), /sine)
            inten[i] = mean( 0.5d * ( (0.5d + c)^2.d + (0.5d + s)^2.d ) )
            if deriv then begin
                dinten[i] = (-1.d/tf) * $
                  mean( (0.5d + c) * cos(!dpi*x^2.d/2.d) + $
                        (0.5d + s) * sin(!dpi*x^2.d/2.d)   ) 
            endif
        endelse
;        print, i, x0, x1, n_cycles, inten[i]
    end

    if n_elements(a) eq 1 then begin
        f = inten
        if deriv then pder = [ [dinten] ]
    endif else begin
        f = a[1] * inten + a[2]
        if deriv then begin
            pder = [ [dinten], [inten], [replicate(1.d,n)] ]
        endif
    endelse

end