#########################################################
#
# Description:
# 	This script performs a variety of simple
#       far field calculations. 
#
#########################################################


# run simulation
clear;
closeall;
run; 

#########################################################
#       Perform the far field projection and image result
E2 = farfield3d("z2"); # this returns |E|^2 in the far field
ux = farfieldux("z2");
uy = farfielduy("z2");

image(ux,uy,E2,"","","|E|^2 at 1 m","polar");


#########################################################
# 	Integrate power in the far field,
#       normalized to the source power 

# pause before running next portion of script
?"Press SPACE BAR to continue";
pause(100); # 

# choose the half angle over which we will integrate
half_angle = 30; #in degrees

# choose the central angle of the cone
cone_center_theta = 0;
cone_center_phi = 0;

# Display the half angle
?"The half angle is: " + num2str(half_angle) + " degrees" +
    " at (theta,phi)=("+num2str(cone_center_theta)+","+num2str(cone_center_phi)+")";


# collect the far field data
E2 = farfield3d("z2"); # this returns |E|^2 in the far field
ux = farfieldux("z2");
uy = farfielduy("z2");


# Method1
temp_cone = farfield3dintegrate(E2,ux,uy,half_angle,cone_center_theta,cone_center_phi);
temp_far  = farfield3dintegrate(E2,ux,uy);
# get fraction of power in cone, then normalize with monitor transmission
T=temp_cone/temp_far*transmission("z2");
?"  The normalized transmission by method 1 is: " + num2str(T*100) + " %";


# Method2
# calculate and integrate the Poynting vector
temp2=0.5*sqrt(eps0/mu0)*farfield3dintegrate(E2,ux,uy,half_angle,cone_center_theta,cone_center_phi);
# apply source normalization 
T=temp2/sourcepower(getdata("z2","f"));

 ?"  The normalized transmission by method 2 is: " + num2str(T*100) + " %";
 
 
#########################################################
# 	Calculate the complex electric field 
#	components in the far field, in cartesian and 
#	spherical coordinate systems

# pause before running next portion of script
?"Press SPACE BAR to continue";
pause(100); # 

# project in cartesian coordinate system
E = farfieldvector3d("z2");
Ex = pinch(E,3,1);
Ey = pinch(E,3,2);
Ez = pinch(E,3,3);
ux = farfieldux("z2");
uy = farfielduy("z2");
image(ux,uy,abs(Ex),"","","|Ex|","polar");
image(ux,uy,abs(Ey),"","","|Ey|","polar");
image(ux,uy,abs(Ez),"","","|Ez|","polar");

# project in spherical (polar) coordinate system
E = farfieldpolar3d("z2",1,201,201);
Er = pinch(E,3,1);
Etheta = pinch(E,3,2);
Ephi = pinch(E,3,3);
ux = farfieldux("z2",1,201,201);
uy = farfielduy("z2",1,201,201);
image(ux,uy,abs(Er),"","","|Er|","polar");
image(ux,uy,abs(Etheta),"","","|Etheta|","polar");
image(ux,uy,abs(Ephi),"","","|Ephi|","polar");



#########################################################
# 	Create line plots of far field data.

# pause before running next portion of script
?"Press SPACE BAR to continue";
pause(100); # 

m="z2";
res = 201;
E2 = farfield3d(m,1,res,res);
ux = farfieldux(m,1,res,res);
uy = farfielduy(m,1,res,res);


# Figure: standard image plot
image(ux,uy,E2,"","","|E|^2","polar");


# Figure: line plot of E2 vs angle at uy=0 (phi=0).  
theta=linspace(-90,90,200);
phi=0;
plot(theta, farfieldspherical(E2,ux,uy,theta,phi),"Theta","E^2 far","Phi="+num2str(phi));


# Figure: plot E2 vs theta at phi=15
# define desired plot
phi=15;
theta=linspace(-90,90,200);
plot(theta,farfieldspherical(E2,ux,uy,theta,phi),"theta","E^2 far","Phi="+num2str(phi));


# Figure: plot E2 vs phi at theta=30
phi=linspace(-180,180,200);
theta=30;
plot(phi,farfieldspherical(E2,ux,uy,theta,phi),"Phi","E^2 far","theta="+num2str(theta));




###########################################
# Calculate the far field disribution 
# at 1m and 10um with the farfield3d and 
# farfieldexact3d functions.  These functions will 
# return the same result in the far field.  
# In the intermediate field, only the 
# farfieldexact3d function is valid.  


# pause before running next portion of script
?"Press SPACE BAR to continue";
pause(100); # 

# set monitor name and resolution
m="z2";
res=101;


# calculate field disribution at 1m.

# calculate various position vectors
ux=farfieldux(m,1,res,1); 
uy=farfielduy(m,1,res,1);
uz=sqrt(1-ux^2-uy^2);
r=1; 
theta=acos(uz)*abs(ux)/ux*180/pi;
if (round(res/2)!=res/2){ # res=odd number, avoid abs(ux)/ux
theta(round(res/2))=0;} # to make sure theta at centre is zero
x=r*sin(theta*pi/180);
y=x*0;
z=r*cos(theta*pi/180);

# Standard far field projection function
E2far1=farfield3d(m,1,res,1);

# Far field exact function
E2exact1=matrix(res);
for (i=1:res) {
  Efar_temp=farfieldexact3d(m,x(i),y(i),z(i));
  E2exact1(i)=sum(abs(Efar_temp)^2);
  if ( (20*round(i/20))==(i) ) { ?"Projection 1: "+num2str(i/res*100)+"%"; }
}

# plot results
plot(theta,E2far1,E2exact1,"theta","|E|^2","1m radius");
legend("farfield3d","farfieldexact3d");



# calculate at distance of 10um
r=10e-6;
x=r*sin(theta*pi/180);
y=x*0;
z=r*cos(theta*pi/180);

# Standard far field projection function
# There is no need to recalculate the projection. 
# Simply scale the fields by 1/r^2
E2far2=E2far1/r^2;

# Far field exact function
E2exact2=matrix(res);
for (i=1:res) {
  Efar_temp=farfieldexact3d(m,x(i),y(i),z(i));
  E2exact2(i)=sum(abs(Efar_temp)^2);
  if ( (20*round(i/20))==(i) ) { ?"Projection 2: "+num2str(i/res*100)+"%"; }
}

plot(theta,E2far2,E2exact2,"theta","|E|^2","10um radius");
legend("farfield3d scaled","farfieldexact3d");


# calculate at a flat plane at x=-10:10,y=0,z=10
x=linspace(-20e-6,20e-6,101);
y=0;
z=10e-6;
E2exact3=farfieldexact3d(m,x,y,z);
E2exact3=abs(E2exact3(1:101,1,1,1)^2) + 
         abs(E2exact3(1:101,1,1,2)^2)^2 + 
         abs(E2exact3(1:101,1,1,3)^2)^2 ;

plot(x*1e6,E2exact3,"x (um)","|E|^2","line at y=0, z=10um");
legend("farfieldexact3d");