%PEC on xz Plane
% test1
% sheet edge parallel to z
% syms a b
% syms fGHz 
% syms lamda k
% syms theta phi 
% syms u v
K=1e3;

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

% 入射波
f0_GHz=0.3;
polAngle_du=0;%电场和入射面夹角
e0t=cosd(polAngle_du); e0p=sind(polAngle_du);
lamda_m=300./f0_GHz./K;
lamda_mm=300./f0_GHz;
k=2.*pi./lamda_m;
theta_iGrid=pi/2;% pi/6入射方向 俯仰角
phi_iGrid_du=linspace(0,180,361);
phi_iGrid=deg2rad(phi_iGrid_du);% 入射方向 方位角
rD0=1e9;

u=sin(theta_iGrid).*cos(phi_iGrid)+sin(theta_iGrid).*cos(phi_iGrid);
v=sin(theta_iGrid).*sin(phi_iGrid)+sin(theta_iGrid).*sin(phi_iGrid);
w=cos(theta_iGrid)+cos(theta_iGrid);

% % Plane xoy E-theta
% 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
Er=-1j.*k.*a.*b./(2.*pi.*rD0).*exp(-1j.*k.*rD0)...
    .*(e0p.*(sin(theta_iGrid).*cos(theta_iGrid).*cos(phi_iGrid)...
    -cos(theta_iGrid).*cos(phi_iGrid).*sin(theta_iGrid))...
    +e0t.*sin(theta_iGrid).*sin(phi_iGrid))...
    .*sinc(k.*u.*a./2).*sinc(k.*w.*b./2);

ErMag=abs(Er); 
ErPhase=angle(Er); ErPhase_du=rad2deg(ErPhase);
ErReal=real(Er);

rcsRTotal=4.*pi.*rD0.^2.*(ErMag).^2;
rcsRTotaldB=10.*log10(rcsRTotal);
% rcsRTotaldB(rcsRTotaldB<-50)=-50;
% figure;
% polarplot(thetaVec,rcsRTotaldB);
% rlim([min(rcsRTotaldB),max(rcsRTotaldB)])
figure;hold on;
plot(phi_iGrid_du,rcsRTotaldB);grid on;
axis tight
figure;
polarplot(phi_iGrid,rcsRTotaldB);grid on;
rlim([min(rcsRTotaldB),max(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;

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

EdGrid=WedgeDiffractionFarField_mono(boolPlot,Ledge,alphaWedge_du,alphaMirrorWedge_du,...
    bdFlag,rEdgePos,phiEdgePos,f0_GHz,rD0,phi_iGrid,theta_iGrid);

% 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_mono(boolPlot,Ledge,alphaWedge_du,alphaMirrorWedge_du,...
    bdFlag,rEdgePos,phiEdgePos,f0_GHz,rD0,phi_iGrid,theta_iGrid);

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

rcsTotal=4.*pi.*rD0.^2.*(ErdMag).^2;
rcsTotaldB=10.*log10(rcsTotal);
figure;
plot(phi_iGrid_du,rcsTotaldB);grid on;
total=[phi_iGrid_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])