001 function [Cb0,Cb1,a0,a1]=rqlf_asympt(u,g,b,S22,S21)
002
003
004
005
006
007
008
009
010
011
012 g=g(:);
013 b=b(:);
014
015 if length(g)~=length(b)
016 error('The lengths of b and g should be the same.')
017 end
018
019
020
021 n=size(S22,1);
022 xtp=mindist(g,b,u,20);
023 xJ=xtp(:);
024 b0=sqrt(sum(xJ.^2));
025 yJ=xJ/b0;
026 PJ=eye(n)-yJ*yJ';
027 dg=-b0*(b+2*(g.*xJ));
028 GJ=-2*b0^2*diag(g)/sqrt(sum(dg.^2));
029 V=S22+S21*S21;
030 sJ2=yJ'*V*yJ;
031 sJ=sqrt(sJ2);
032 a0=sJ*sqrt(1-yJ'*S21*(eye(n)+GJ)*S21*yJ/sJ2)/sqrt(det(eye(n)+ ...
033 PJ*GJ*PJ));
034 aJ=-1/sJ*PJ*(eye(n)+GJ)*S21*yJ;
035 WJ1=inv(eye(n)+PJ*GJ*PJ+aJ*aJ');
036 WJ2=inv(eye(n)+PJ*GJ*PJ);
037 vJ1=PJ*S21*PJ*GJ*yJ;
038 vJ2=PJ*GJ*yJ;
039 vJ3=PJ*GJ*S21*yJ;
040 vJ4=PJ*GJ*PJ*V*yJ;
041 KJ1=yJ'*S21*GJ*yJ;
042 KJ2=yJ'*GJ*yJ;
043 QJ1=PJ*GJ*PJ;
044 QJ2=PJ*S21*PJ*GJ*PJ;
045 QJ3=QJ2+.5*KJ1*QJ1-vJ2*vJ3';
046 QJ4=-(yJ'*V*PJ*GJ*yJ)*PJ*GJ*PJ-2*PJ*GJ*PJ*V*yJ*yJ'*GJ*PJ+PJ*GJ* ...
047 PJ*V*PJ*GJ*PJ;
048 a11=1/2*holmquist1(WJ1,QJ4)/sJ2-1/2*holmquist1(WJ1,vJ4*vJ4')/sJ2^2+ ...
049 1/2*holmquist2(WJ1,QJ3,QJ3)/sJ2+holmquist2(WJ1,aJ*vJ4',QJ3)/ ...
050 sJ^3+1/2*holmquist2(WJ1,aJ*aJ',vJ4*vJ4')/sJ2^2+1/8* ...
051 holmquist3(WJ1,vJ2*vJ2',QJ1,QJ1)-1/2*holmquist2(WJ1,vJ2*vJ2',QJ1)-1/8*(KJ2+1)*holmquist2(WJ1,QJ1,QJ1)+...
052 1/2*holmquist2(WJ1,vJ2*vJ4',QJ1)/sJ2;
053 a12=(1/2*holmberg1(WJ2,aJ,vJ1,QJ1)+holmberg1(WJ2,aJ,vJ2,QJ2)-1/2*holmberg1(WJ2,aJ,vJ3,QJ1)+...
054 KJ1*holmberg1(WJ2,aJ,vJ2,QJ1)-holmberg1(WJ2,aJ,vJ3,vJ2*vJ2')-1/2*KJ2*holmberg1(WJ2,aJ,vJ3,QJ1))/sJ-...
055 1/2*holmberg2(WJ2,aJ,vJ2,QJ3,QJ1)/sJ-1/2*holmberg2(WJ2,aJ,aJ,vJ2*vJ2',QJ1)-...
056 1/8*(KJ2+1)*holmberg2(WJ2,aJ,aJ,QJ1,QJ1)+1/8*holmberg3(WJ2,aJ,aJ,vJ2*vJ2',QJ1,QJ1);
057 [det(WJ1) det(WJ2) sJ a11 a12];
058 a1=sJ*(sqrt(det(WJ1)*a11)+sqrt(det(WJ2))*a12);
059 Cb0=a0*exp(-b0^2/2)/pi;
060 Cb1=a1*exp(-b0^2/2)/pi/b0^2;
061
062
063
064
065
066
067
068