001 function bvn = bvnormcdf( b1, b2, r )
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045 error(nargchk(3,3,nargin))
046 [errorcode,h,k,r] = comnsize(-b1,-b2,r);
047 if (errorcode > 0)
048 error ('b1, b2 and r must be of common size or scalar');
049 end
050
051 bvn = zeros(size(h));
052
053 zero = 0.d0;
054 one = 1.d0;
055 two = 2.d0;
056 twopi = 6.283185307179586d0;
057
058 w6 = [ 0.1713244923791705D+00, 0.3607615730481384D+00, 0.4679139345726904D+00];
059 x6 = [-0.9324695142031522D+00,-0.6612093864662647D+00,-0.2386191860831970D+00];
060
061 w12 = [ 0.4717533638651177d-01, 0.1069393259953183d+00, 0.1600783285433464d+00,...
062 0.2031674267230659d+00, 0.2334925365383547d+00, 0.2491470458134029d+00];
063 x12=[ -0.9815606342467191d+00,-0.9041172563704750d+00,-0.7699026741943050d+00, ...
064 -0.5873179542866171d+00,-0.3678314989981802d+00,-0.1252334085114692d+00];
065
066 w20 = [ 0.1761400713915212d-01, 0.4060142980038694d-01, ...
067 0.6267204833410906d-01, 0.8327674157670475d-01, ...
068 0.1019301198172404d+00, 0.1181945319615184d+00,...
069 0.1316886384491766d+00, 0.1420961093183821d+00,...
070 0.1491729864726037d+00, 0.1527533871307259d+00];
071 x20=[ -0.9931285991850949d+00, -0.9639719272779138d+00, ...
072 -0.9122344282513259d+00, -0.8391169718222188d+00, ...
073 -0.7463319064601508d+00,-0.6360536807265150d+00,...
074 -0.5108670019508271d+00,-0.3737060887154196d+00, ...
075 -0.2277858511416451d+00, -0.7652652113349733d-01];
076
077 hk = h.*k;
078
079 k0 = find(abs(r) < 0.925d0 );
080 if (k0)
081 hs = ( h(k0).^2 + k(k0).^2 )/two;
082 asr = asin(r(k0));
083 k1 = find(r(k0)>0.75);
084 if any(k1)
085 k01 = k0(k1);
086 for i = 1:10
087 for is = -1:2:1,
088 sn = sin( asr(k1).*(is.*x20(i) +1 )/2 );
089 bvn(k01) = bvn(k01) + w20(i)*exp( ( sn.*hk(k01) - hs(k1) )./(1 - sn.*sn ) );
090 end
091 end
092 end
093 k1 = find(0.3 <= r(k0) & r(k0)<0.75);
094 if any(k1)
095 k01 = k0(k1);
096 for i = 1:6
097 for is = -1:2:1,
098 sn = sin( asr(k1).*(is.*x12(i) +1 )/2 );
099 bvn(k01) = bvn(k01) + w12(i)*exp( ( sn.*hk(k01) - hs(k1) )./(1 - sn.*sn ) );
100 end
101 end
102 end
103 k1 = find( r(k0)<0.3);
104 if any(k1)
105 k01 = k0(k1);
106 for i = 1:3
107 for is = -1:2:1,
108 sn = sin( asr(k1).*(is.*x6(i) +1 )/2 );
109 bvn(k01) = bvn(k01) + w6(i)*exp( ( sn.*hk(k01) - hs(k1) )./(1 - sn.*sn ) );
110 end
111 end
112 end
113 bvn(k0) = bvn(k0).*asr/( two*twopi ) + fi(-h(k0)).*fi(-k(k0));
114 end
115
116 k1 = find(0.925<=abs(r) & abs(r)<=1 );
117 if any(k1)
118 k2 = find(r(k1) < 0);
119 if any(k2 )
120 k12 = k1(k2);
121 k(k12) = -k(k12);
122 hk(k12) = -hk(k12);
123 end
124 k3 = find( abs(r(k1)) < 1);
125 if (k3)
126 k13 = k1(k3);
127 as = (1 - r(k13) ).*(1 + r(k13) );
128 a = sqrt(as);
129 b = abs( h(k13) - k(k13) );
130 bs = b.*b;
131 c = ( 4.d0 - hk(k13) )/8.d0;
132 d = ( 12.d0 - hk(k13) )/16.d0;
133 asr = -( bs./as + hk(k13) )/2.d0;
134 k4 = find(asr > -100.d0);
135 if (k4)
136 bvn(k13(k4)) = a(k4).*exp(asr(k4)).*(1 - c(k4).* ....
137 ( bs(k4) - as(k4)).*(1 - d(k4).*bs(k4)/5 )/3 ...
138 + c(k4).*d(k4).*as(k4).^2/5 );
139 end
140 k5 = find(hk(k13) < 100.d0);
141 if ( k5 )
142
143 k135 = k13(k5);
144 bvn(k135) = bvn(k135) - exp( -hk(k135)/2 ).*sqrt(twopi).*fi(-b(k5)./a(k5)).*b(k5)...
145 .*(1- c(k5).*bs(k5).*(1 - d(k5).*bs(k5)/5 )/3 );
146 end
147 a = a/two;
148 for i = 1:10
149 for is = -1:2:1,
150 xs = ( a.*( is*x20(i) + 1 ) ).^2;
151 rs = sqrt(1 - xs );
152 asr = -( bs./xs + hk(k13) )/2;
153 k6 = find( asr > -100.d0 ) ;
154 if any(k6)
155 k136 = k13(k6);
156 bvn(k136) = bvn(k136) + a(k6).*w20(i).*exp( asr(k6) )...
157 .*( exp(-hk(k136).*(1-rs(k6))./(2*(1+rs(k6))))./rs(k6)...
158 - (1 + c(k6).*xs(k6).*(1 + d(k6).*xs(k6) ) ) );
159 end
160 end
161 end
162 bvn(k3) = -bvn(k3)/twopi;
163
164 end
165 k7 = find(r(k1)>0);
166 if any(k7 ),
167 k17 = k1(k7);
168 bvn(k17) = bvn(k17) + fi( -max( h(k17), k(k17)) );
169 end
170 k8 = find(r(k1)<0);
171 if any(k8),
172 k18 = k1(k8);
173 bvn(k18) = -bvn(k18) + max(0, fi(-h(k18)) - fi(-k(k18)) );
174 end
175 end
176
177 k9 = find(abs(r)>1);
178 if any(k9)
179 bvn(k9) = NaN;
180 end
181 val = bvn;
182 return
183
184 function F = fi(x)
185 F = 0.5.*(erfc((-x)./sqrt(2)));
186
187