Add files via upload

PhC2D_hex_PWE_f.m

@@ -0,0 +1,86 @@
+function[E,f0]=PhC2D_hex_PWE_f(Xhex,Yhex,Gxhex,Gyhex,k,HHH,nmodes,TE,TM);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+c=2.99792458e8;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%% Interpolation on a grid that have 2^N points %%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nx=length(Xhex(1,:));
+Ny=length(Yhex(:,1));
+
+NGx=length(Gxhex(1,:));
+NGy=length(Gyhex(:,1));
+NG=NGx*NGy;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%% Building Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+GXX=Gxhex(:);
+GYY=Gyhex(:);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if TE==1
+    GkXX = ( GXX + k(1) )*( GXX + k(1) )'; % Gk(i,j) = (G(i) + k)*(G(j) + k)
+    GkYY = ( GYY + k(2) )*( GYY + k(2) )';
+    Gk=GkXX+GkYY;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if TM==1
+    Gk1 = sqrt( ( GXX + k(1) ) .^2 + ( GYY + k(2) ) .^2 )  ;
+    Gk2 = sqrt( ( GXX + k(1) )'.^2 + ( GYY + k(2) )'.^2 );
+    Gk=Gk1*Gk2;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+HH=reshape(Gk,[NGy,NGx,NGy,NGx]);
+
+H=HH.*HHH;
+H=reshape(H,NG,NG);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% Solving Hamiltonian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[psik,k0] = eig(H);   %% eigen values are ordered
+
+f0 = sqrt(diag(k0)) ; % actually it is w0
+lambda= 2*pi ./ sqrt(diag(k0)) * 1e6  ;
+
+f0=f0(1:nmodes);
+psik = psik(:,1:nmodes);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%% Reverse Hexagonal Fourier Transform %%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for j=1:nmodes
+
+PSI = reshape(psik(:,j),[NGy,NGx]);
+Ghex=PSI;
+
+for m=1:Nx
+for l=1:Ny
+    wwhex = exp(1i*(   (Xhex(l,m)-Xhex(1))*(Nx-1)/Nx *Gxhex  +  (Yhex(l,m)-Yhex(1))*(Ny-1)/Ny *Gyhex   )  );
+    Fhex(l,m) = (1/(Nx*Ny))*sum(sum(Ghex.*wwhex));
+end
+end
+
+    E(:,:,j)= Fhex /max(Fhex(:)); 
+end
+
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PhC2D_sq_PWE_f.m

@@ -0,0 +1,97 @@
+function[E,f0]=PhC2D_sq_PWE_f(x,y,Gx,Gy,k,HHH,nmodes,TE,TM);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+c=2.99792458e8;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%% Interpolation on a grid that have 2^N points %%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nx=length(x);
+Ny=length(y);
+
+NGx=length(Gx);
+NGy=length(Gy);
+NG=NGx*NGy;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%% Building Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[GXX,GYY]=meshgrid(Gx,Gy);
+GXX=GXX(:);
+GYY=GYY(:);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if TE==1
+    GkXX = ( GXX + k(1) )*( GXX + k(1) )'; % Gk(i,j) = (G(i) + k)*(G(j) + k)
+    GkYY = ( GYY + k(2) )*( GYY + k(2) )';
+    Gk=GkXX+GkYY;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if TM==1
+    Gk1 = sqrt( ( GXX + k(1) ) .^2 + ( GYY + k(2) ) .^2 )  ;
+    Gk2 = sqrt( ( GXX + k(1) )'.^2 + ( GYY + k(2) )'.^2 );
+    Gk=Gk1*Gk2;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+HH=reshape(Gk,[NGy,NGx,NGy,NGx]);
+
+H=HH.*HHH;
+H=reshape(H,NG,NG);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% Solving Hamiltonian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[psik,k0] = eig(H);   %% eigen values are ordered
+
+f0 = sqrt(diag(k0)) ; % actually it is w0
+lambda= 2*pi ./ sqrt(diag(k0)) * 1e6  ;
+
+f0=f0(1:nmodes);
+psik = psik(:,1:nmodes);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%% Transforming & Scaling the waves functions %%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for j=1:nmodes
+    PSI = reshape(psik(:,j),[NGy,NGx]);
+    PSI = invFFT2D(PSI,Ny,Nx);
+    E(:,:,j) = PSI / max(PSI(:));
+end
+
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [Vxy] = invFFT2D(Vk2D,Ny,Nx)
+
+Nkx=length(Vk2D(1,:));
+Nky=length(Vk2D(:,1));
+
+Nx1=Nx/2-floor(Nkx/2);
+Nx2=Nx/2+ceil(Nkx/2);
+Ny1=Ny/2-floor(Nky/2);
+Ny2=Ny/2+ceil(Nky/2);
+
+Vk2D00=zeros(Ny,Nx);
+Vk2D00( Ny1+1:Ny2 , Nx1+1:Nx2)=Vk2D;
+Vxy=ifft2(ifftshift(Vk2D00));
+
+end