0001 function xyredfac = redfac_dke_jd(waveparam,xyD0_rf_mod,betath_ref,limit_dql);
0002
0003
0004 Nparmin = -10;
0005 Nparmax = 10;
0006 nz = 1000;
0007
0008 nit_max = 5000;
0009
0010 tol = 1e-10;
0011
0012 [nr,nb] = size(xyD0_rf_mod);
0013
0014 xyredfac = ones(nr,nb);
0015
0016 for ir = 1:nr,
0017
0018 ydNpar = waveparam(ir,:,8);
0019 yNpar0 = waveparam(ir,:,7);
0020 yepolz = waveparam(ir,:,12);
0021 yD0_rf_mod = xyD0_rf_mod(ir,:);
0022
0023 yalphab = betath_ref*abs(yepolz).^2.*abs(yNpar0).*yD0_rf_mod./ydNpar/sqrt(pi);
0024 zNpar = linspace(Nparmin,Nparmax,nz);
0025
0026 zyalphab = repmat(yalphab,[nz,1]);
0027 zyNpar = repmat(zNpar',[1,nb]);
0028 zydNpar = repmat(ydNpar,[nz,1]);
0029 zyNpar0 = repmat(yNpar0,[nz,1]);
0030
0031 zyDparparb = zyalphab.*exp(-(zyNpar - zyNpar0).^2./zydNpar.^2);
0032
0033 yredfac = ones(1,nb);
0034
0035 it_rfac = 0;
0036
0037 while it_rfac < nit_max,
0038
0039 it_rfac = it_rfac + 1;
0040
0041 zyredfac = repmat(yredfac,[nz,1]);
0042 zyDparparb_mod = zyDparparb./zyredfac;
0043 zDparpar = sum(zyDparparb_mod,2);
0044 Dmax = max(zDparpar);
0045 if Dmax - limit_dql < tol,
0046 break
0047 end
0048 izmax = min(find(Dmax == zDparpar));
0049 Dmaxmax = max(zyDparparb_mod(izmax,:));
0050 iymax = find(abs(zyDparparb_mod(izmax,:) - Dmaxmax) < tol);
0051 iymin = find(abs(zyDparparb_mod(izmax,:) - Dmaxmax) >= tol);
0052 if ~isempty(iymin);
0053 Dmaxmin = max(zyDparparb_mod(izmax,iymin));
0054 if sum(zyDparparb_mod(izmax,iymin),2) + Dmaxmin*length(iymax) > limit_dql,
0055 Dmaxnew = Dmaxmin;
0056 else
0057 Dmaxnew = (limit_dql - sum(zyDparparb_mod(izmax,iymin),2))/length(iymax);
0058 end
0059 else
0060 Dmaxnew = limit_dql/nb;
0061 end
0062 yredfac(iymax) = yredfac(iymax).*zyDparparb_mod(izmax,iymax)/Dmaxnew;
0063 end
0064
0065 if it_rfac == nit_max,
0066 error('the calculation of xyredfac failed')
0067 end
0068
0069 xyredfac(ir,:) = yredfac;
0070 end
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086