# Written in Julia by André Ferrari
# January 2015
#
# Translation of matlab files bvn.m and bvnu.m authored by:
#
# Alan Genz
# Department of Mathematics
# Washington State University
# Pullman, WA 99164-3113
# Email : alangenz@wsu.edu
#
# Original source available from
# http://www.math.wsu.edu/faculty/genz/software/matlab/bvn.m
#
# Copyright (C) 2013, Alan Genz,  All rights reserved.               
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided the following conditions are met:
#   1. Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#   2. Redistributions in binary form must reproduce the above copyright
#      notice, this list of conditions and the following disclaimer in 
#      the documentation and/or other materials provided with the 
#      distribution.
#   3. The contributor name(s) may not be used to endorse or promote 
#      products derived from this software without specific prior 
#      written permission.
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS 
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR 
# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#

function bvn(xl::Number,xu::Number,yl::Number,yu::Number,r::Number)
    p = bvnu(xl,yl,r) - bvnu(xu,yl,r) - bvnu(xl,yu,r) + bvnu(xu,yu,r)
    p = max(0.0, min(p,1.0))
end

phid(z::Number) = erfc(-z/sqrt(2.0))/2.0
    
function bvnu(dh::Number,dk::Number,r::Number)
    if dh ==  Inf || dk ==  Inf
        p = 0.0
    elseif dh == -Inf 
        if dk == -Inf 
            p = 1.0
        else 
            p = phid(-dk)
        end
    elseif dk == -Inf
        p = phid(-dh)
    elseif r == 0.0
        p = phid(-dh)*phid(-dk)  
    else 
        tp = 2.0*pi
        h = dh
        k = dk
        hk = h*k
        bvn = 0.0
        if abs(r) < 0.3      
            w = [0.1713244923791705, 0.3607615730481384, 0.4679139345726904]
            x = [0.9324695142031522, 0.6612093864662647, 0.2386191860831970]
        elseif abs(r) < 0.75 
            w = [.04717533638651177, 0.1069393259953183, 0.1600783285433464]
            w = [w, 0.2031674267230659, 0.2334925365383547, 0.2491470458134029]
            x = [0.9815606342467191, 0.9041172563704750, 0.7699026741943050]
            x = [x, 0.5873179542866171, 0.3678314989981802, 0.1252334085114692]
        else                
            w = [.01761400713915212, .04060142980038694, .06267204833410906]
            w = [w, .08327674157670475, 0.1019301198172404, 0.1181945319615184]
            w = [w, 0.1316886384491766, 0.1420961093183821, 0.1491729864726037]
            w = [w, 0.1527533871307259]
            x = [0.9931285991850949, 0.9639719272779138, 0.9122344282513259]
            x = [x, 0.8391169718222188, 0.7463319064601508, 0.6360536807265150]
            x = [x, 0.5108670019508271, 0.3737060887154196, 0.2277858511416451]
            x = [x, 0.07652652113349733]
        end 
        w = [w, w]
        x = [1.0-x, 1.0+x]
        
        if abs(r) < 0.925             
            hs = ( h*h + k*k )/2
            asr = asin(r)/2.0  
            sn = sin(asr*x)
            bvn = dot(exp((sn*hk-hs)./(1-sn.^2)),w)
            bvn = bvn*asr/tp + phid(-h)*phid(-k)            
        else             
            if r < 0.0
                k = -k
                hk = -hk
            end
            if abs(r) < 1.0
                as = 1.0-r^2
                a = sqrt(as)
                bs = (h-k)^2
                asr = -( bs/as + hk )/2.0
                c = (4.0-hk)/8.0 
                d = (12.0-hk)/80.0
                if asr > -100.0
                    bvn = a*exp(asr)*(1.0-c*(bs-as)*(1-d*bs)/3.0+c*d*as^2)
                end
                if hk  > -100.0
                    b = sqrt(bs)
                    sp = sqrt(tp)*phid(-b/a)
                    bvn = bvn - exp(-hk/2.0)*sp*b*( 1.0 - c*bs*(1-d*bs)/3.0 )
                end   
                a = a/2.0
                xs = (a*x).^2
                asr = -( bs./xs + hk )/2.0
                ix = find( asr .> -100.0 )
                xs = xs[ix]
                sp = ( 1.0 + c*xs.*(1.0+5.0*d*xs))
                rs = sqrt(1.0-xs)
                ep = exp( -(hk/2.0)*xs./(1+rs).^2 )./rs
                bvn = (a*dot( exp(asr[ix]).*(sp-ep),w[ix]) - bvn)/tp                
            end 

            if r > 0.0
                bvn =  bvn + phid(-max(h,k))
            elseif h >= k
                bvn = -bvn
            else 
                if h < 0.0
                    L = phid(k)-phid(h)
                else
                    L = phid(-h)-phid(-k)
                end
                bvn =  L - bvn
            end
        end
        p = max(0.,min(1.,bvn))
    end
    
    return p
    
end