function F=pwd(u,rWD,RD,Sf1,Sf2,C1D,C2D,wuf11,wuf21,lanmf11,lanmf21,n11D,n21D,n12D,n22D,M112,M212,tafa11,tafa21,Pi1D,Pi2D)
f1u=(1-wuf11)/(lanmf11*n11D+(1-wuf11)*u);
f2u=(1-wuf21)/(lanmf21*n21D+(1-wuf21)*u);
a11=rWD*sqrt(f1u)*besseli(1,rWD*sqrt(f1u));
a12=-rWD*sqrt(f1u)*besselk(1,rWD*sqrt(f1u));
a21=besseli(0,RD*sqrt(f1u));
a22=besselk(0,RD*sqrt(f1u));
a23=-besselk(0,RD*sqrt(1/n12D));
a31=M112*sqrt(f1u)*besseli(1,RD*sqrt(f1u));
a32=-M112*sqrt(f1u)*besselk(1,RD*sqrt(f1u));
a33=sqrt(1/n12D)*besselk(1,RD*sqrt(1/n12D));
b11=rWD*sqrt(f2u)*besseli(1,rWD*sqrt(f2u));
b12=-rWD*sqrt(f2u)*besselk(1,rWD*sqrt(f2u));
b21=besseli(0,RD*sqrt(f2u));
b22=besselk(0,RD*sqrt(f2u));
b23=-besselk(0,RD*sqrt(1/n22D));
b31=M212*sqrt(f2u)*besseli(1,RD*sqrt(f2u));
b32=-M212*sqrt(f2u)*besselk(1,RD*sqrt(f2u));
b33=sqrt(1/n22D)*besselk(1,RD*sqrt(1/n22D));
A11=-(a22*a33-a23*a32)/(tafa11*u*(a11*a22*a33-a11*a23*a32-a12*a21*a33+a12*a23*a31));
B11=(a21*a33-a23*a31)/(tafa11*u*(a11*a22*a33-a11*a23*a32-a12*a21*a33+a12*a23*a31));
A21=-(b22*b33-b23*b32)/(tafa21*u*(b11*b22*b33-b11*b23*b32-b12*b21*b33+b12*b23*b31));
B21=(b21*b33-b23*b31)/(tafa21*u*(b11*b22*b33-b11*b23*b32-b12*b21*b33+b12*b23*b31));
X1=@(aerfa)(A11*besseli(0,sqrt(f1u*(0.732-aerfa).^2))+B11*besselk(0,sqrt(f1u*(0.732-aerfa).^2)));
X2=@(aerfa)(A21*besseli(0,sqrt(f2u*(0.732-aerfa).^2))+B21*besselk(0,sqrt(f2u*(0.732-aerfa).^2)));
Y1=quadgk(X1,-1,1,'RelTol',1e-8,'AbsTol',1e-12);
Y2=quadgk(X2,-1,1,'RelTol',1e-8,'AbsTol',1e-12);
puw1D=Y1/2+Sf1/u+Pi1D/u;
puw2D=Y2/2+Sf2/u+Pi2D/u;
pwD=((1/(puw1D*u.^2-u*Pi1D)+C1D)*Pi1D+(1/(puw2D*u.^2-u*Pi2D)+C2D)*Pi2D+1/u)/(u*(C1D+C1D)+(1/(u*puw1D-Pi1D)+1/(u*puw2D-Pi2D)));
F=pwD;
end