%PEC on xz Plane
% sheet edge parallel to z
% syms a b
% syms fGHz 
% syms lamda k
% syms theta phi 
% syms u v
K=1e3;du2rad=pi./180;eps=1e-6;

a=5;%m Dx 边长
b=5;%m Dy

% 入射波
f0_GHz=0.3;
polAngle_du=90;%电场和入射面夹角
e0t=cosd(polAngle_du); e0p=sind(polAngle_du);
lamda_m=300./f0_GHz./K;
lamda_mm=300./f0_GHz;
k=2.*pi./lamda_m;
thetaInc_i0=pi/2;% pi/6入射方向 俯仰角
phiInc_i0=pi/4;% 入射方向 方位角

% 观察角度
nr=501; nTheta=361;
rDGrid=1000; thetaD=pi/2;%linspace(-pi,pi,nTheta)
phiVec_du=linspace(0,360,nTheta); phiVec=deg2rad(phiVec_du);
phiDGrid=phiVec; %phiDGrid=0;

xGrid=rDGrid.*cos(phiDGrid); yGrid=rDGrid.*sin(phiDGrid);
xGrid_mm=xGrid.*1000; yGrid_mm=yGrid.*1000;

u=sin(thetaD).*cos(phiDGrid)+sin(thetaInc_i0).*cos(phiInc_i0);
v=sin(thetaD).*sin(phiDGrid)+sin(thetaInc_i0).*sin(phiInc_i0);
w=cos(thetaD)+cos(thetaInc_i0);

% % Plane xoy E-theta
% e0t=1; e0p=0;
% Er=-1j.*k.*a.*b./(2.*pi.*rDGrid).*exp(-1j.*k.*rDGrid).*cos(thetaD)...
%     .*(e0t.*cos(phiDGrid-phiInc_i0)...
%     +e0p.*cos(thetaInc_i0).*sin(phiDGrid-phiInc_i0))...
%     .*sinc(k.*a.*u./2).*sinc(k.*b.*v./2);

% Plane xoz E-theta
% e0t=1; e0p=0;
Ert=-1j.*k.*a.*b./(2.*pi.*rDGrid).*exp(-1j.*k.*rDGrid)...
    .*(e0p.*(sin(thetaD).*cos(thetaInc_i0).*cos(phiInc_i0)...
    -cos(thetaD).*cos(phiDGrid).*sin(thetaInc_i0))...
    +e0t.*sin(thetaD).*sin(phiInc_i0))...
    .*sinc(k.*u.*a./2).*sinc(k.*w.*b./2);
Erp=-1j.*k.*a.*b./(2.*pi.*rDGrid).*exp(-1j.*k.*rDGrid)...
    .*e0p.*sin(phiDGrid).*sin(thetaInc_i0)...
    .*sinc(k.*u.*a./2).*sinc(k.*w.*b./2);
Er=(e0t==1 && e0p==0).*Ert+(e0t==0 && e0p==1).*Erp;
ErMag=abs(Er); 
ErPhase=angle(Er); ErPhase_du=rad2deg(ErPhase);
ErReal=real(Er);

rcsRTotal=4.*pi.*rDGrid.^2.*(ErMag).^2;
rcsRTotaldB=10.*log10(rcsRTotal);
% rcsRTotaldB(rcsRTotaldB<-50)=-50;
% figure;
% polarplot(thetaVec,rcsRTotaldB);
% rlim([min(rcsRTotaldB),max(rcsRTotaldB)])
figure;%(1)hold on;
plot(phiVec_du,rcsRTotaldB);grid on;
figure;%hold on;%(2)
% figure(2);hold on;%
polarplot(phiVec,rcsRTotaldB);grid on;
% rlim([min(rcsRTotaldB),max(rcsRTotaldB)])
rlim([-20,max(rcsRTotaldB)])

totalOp=[phiVec_du;rcsRTotaldB]';
% figure;
% polarplot(thetaVec,rcsRTotal);
% figure;
% plot3D(xGrid,yGrid,ErMag,ErMag);
% title('ErMag');axis equal;% axis tight;
% figure;
% plot3D(xGrid,yGrid,ErPhase,ErPhase);
% title('ErPhase');axis equal;% axis tight;
% figure;
% plot3D(xGrid,yGrid,ErReal,ErReal);
% title('ErReal');axis equal;% axis tight;

%%
beta0_du=phiInc_i0-(-90);%入射线和边缘夹角
bdFlag=-1;%边界类型 -1 狄拉克边界 边界值为0; 1 黎曼边界 边界法向导数为0
alphaWedge_du=0;%楔形边缘夹角
alphaMirrorWedge_du=-180;%楔形角的对称位置
Ledge=5;%边缘长度
rEdgePos=a/2; %边缘位置
phiEdgePos=0;
boolPlot=0;%画图

EdGrid=WedgeDiffractionFarField(boolPlot,Ledge,alphaWedge_du,alphaMirrorWedge_du,...
    bdFlag,rEdgePos,phiEdgePos,f0_GHz,rDGrid,phiDGrid,thetaD,phiInc_i0,thetaInc_i0);

% EdMag=abs(Ed); 
% rcsdTotal=4.*pi.*rGrid.^2.*(EdMag).^2;
% rcsdTotaldB=10.*log10(rcsdTotal);
% figure;
% plot(thetaVec_du,rcsdTotaldB);

alphaMirrorWedge_du=0;%楔形角的对称位置
rEdgePos=-a/2; %边缘位置
phiEdgePos=pi;
EdGrid1=WedgeDiffractionFarField(boolPlot,Ledge,alphaWedge_du,alphaMirrorWedge_du,...
    bdFlag,rEdgePos,phiEdgePos,f0_GHz,rDGrid,phiDGrid,thetaD,phiInc_i0,thetaInc_i0);

Erd=Er+EdGrid+EdGrid1;
% Erd=0+EdGrid+EdGrid1;
ErdMag=abs(Erd); 
ErdPhase=angle(Erd); ErdPhase_du=rad2deg(ErdPhase);
ErdReal=real(Erd);

rcsTotal=4.*pi.*rDGrid.^2.*(ErdMag).^2;
rcsTotaldB=10.*log10(rcsTotal);
figure(4);hold on;
plot(phiVec_du,rcsTotaldB);grid on;
total=[phiVec_du;rcsTotaldB]';
% figure;
% polarplot(thetaVec,rcsTotaldB);
% figure;
% polarplot(thetaVec,rcsTotal);
% figure;
% plot3D(xGrid,yGrid,ErdMag,ErdMag);
% title('ErMag');axis equal;% axis tight;
% figure;
% plot3D(xGrid,yGrid,ErdPhase,ErdPhase);
% title('ErPhase');axis equal;% axis tight;
% figure;
% plot3D(xGrid,yGrid,ErdReal,ErdReal);
% title('ErReal');axis equal;% axis tight;

% figure(13)
% caxis([0 20])