# script file: mwp_45deg.lsf
#
# This file modifies mwp.fsp, rotating the source by 45 degrees.  
# A number of images are created, and the tranmission is calculated.
# The simulation volume is increase to make better field images, so
# several periods are visible. If you are only interested in 
# transmission measurements, a smaller simulation volume is better.


# if(1)	- open mwp.fsp
# 	- modify setup
# 	- run simulation
# if(0)	- skip the above steps, go straight to the analysis section
if(1) {
	runtype = 1; # set to 2 for non-graphical, 3 for parallel
	load("mwp");
	save("mwp45.fsp");	#create a new file for these simulations

	switchtolayout;
	W=70e-9;	 # wire width	
	setnamed("::model","grating_width",W);
	setnamed("::model","duty_cycle",0.5);	# 50% duty cycle

	select("FDTD");
	set("y", 0.25e-6);
	set("y span",1.5e-6);
	set("x min bc","Bloch");		#must use Bloch boundary conditions
	set("pml reflection",1e-10);			#large angles require smaller reflection
	set("set based on source angle",1);
	set("simulation time",1000e-15);
	set("mesh accuracy",4);
	select("source1");
	set("angle",45);
	set("y",0.90e-6);			#move source close to edge
	lambda=550e-9; 
	set("center wavelength",lambda);
	set("wavelength span",0);
	set("polarization angle",0);

	setnamed("reflection","y",0.95e-6);	#move old monitors close to the edge
	setnamed("t2","y",0.95e-6);		
	setnamed("transmission","y",-0.45e-6);
	setnamed("t1","y",-0.45e-6);	
	addprofile;
	set("name", "fields");
	set("monitor type","2D z-normal");	#2d
	set("x",0);
	set("y",0);
	set("x span",W*3);
	set("y span",2.5e-6);

	# Run simulation	
	run(runtype);
}

# Analysis section
select("FDTD");
kx=get("kx");

Ex=pinch(getdata("fields","Ex"));
Ey=pinch(getdata("fields","Ey"));
Emag=pinch(getelectric("fields"));
x=getdata("fields","x");
y=getdata("fields","y");

# Expand the data into 5 unit cells, to make the graphs look better.
# First must make 5 copies of the data.  
mat_size=size(Ex);
n1=mat_size(1);
n2=mat_size(2);

Emag2=matrix(n1*5,n2);
Ex2=matrix(n1*5,n2);
Ey2=matrix(n1*5,n2);
x2= matrix(n1*5,1);

for (cnt=0:4) {
  indices_x = (1:n1)+cnt*n1;
  indices_y = (1:n2);

  Emag2(indices_x,indices_y) = Emag;
  Ex2(  indices_x,indices_y) = Ex * exp(1i*cnt*kx*get("x span"));
  Ey2(  indices_x,indices_y) = Ey * exp(1i*cnt*kx*get("x span"));
  x2(   indices_x)           = x + cnt * (x(n1)-x(1));

}

image(x2,y,Emag2, "x","y","|E|^2");
setplot("colorbar min",0); setplot("colorbar max",3);
image(x2,y,real(Ex2),"x","y","real(Ex)");
setplot("colorbar min",-1); setplot("colorbar max",1);
image(x2,y,real(Ey2),"x","y","real(Ey)");
setplot("colorbar min",-1); setplot("colorbar max",1);
image(x2,y,angle(Ex2),"x","y","Ex Phase");
image(x2,y,angle(Ey2),"x","y","Ey Phase");

?"total transmission is "+num2str(abs(transmission("transmission")));	#print out the transmission