##############################################################
# Scriptfile: OLED_hex_analysis
# This file calculates farfield distribution from
# dipoles distribute randomly in the unit cell 
# using the mirror symmetry and the 60 degree rotational symmetry.
# Also, this calculates enhancement of extraction efficiency
#
# In the far_field_change_index analysis group
# E2_far1 is glass
# E2_far2 is air  
############################################################## 

#this resolution can greatly affect the time needed to run this script
nthe=91;
nphi=181; #nphi should be 6*n+1 (n;integer);

save_data=1; # save data to a ldf file

loaddata("num_locs");

#extraction of frequency data
load("OLED_hex_simple0_x");
temp = getresult("far_field_change_index","E2_far");
save;
f = temp.f;
na = length(temp.ux);
nb = length(temp.uy);

E2x=matrix(na,nb,length(f));
E2y=matrix(na,nb,length(f));
E2z=matrix(na,nb,length(f));
E2sp_tot=matrix(nthe,nphi,length(f));

##### calculation of far field radiated from dipoles distributed randomly in the unit cell with PC -start-
##### calculation for regions from 2 to 6 -start- 
E2spA_tot=matrix(nthe,nphi,length(f));

if(ireg_1p) {
for(nn=1:ireg_1p) {

load("OLED_hex_simple_"+"1p"+num2str(nn)+"_x");
?"1p, x, nn="+num2str(nn);
E2_far=getresult("far_field_change_index","E2_far");
save;
ux=E2_far.ux;
uy=E2_far.uy;
E2x=E2_far.E2_far1;

load("OLED_hex_simple_"+"1p"+num2str(nn)+"_y");
?"1p, y, nn="+num2str(nn);
E2_far=getresult("far_field_change_index","E2_far");
save;
E2y=E2_far.E2_far1;

load("OLED_hex_simple_"+"1p"+num2str(nn)+"_z");
?"1p, z, nn="+num2str(nn);
E2_far=getresult("far_field_change_index","E2_far");
save;
E2z=E2_far.E2_far1;

E2=E2x+E2y+E2z; #far field from dipoles located in region 1'. 

#transfomation of coordinate from (ux,uy) to (theta,phi) using farfieldspherical function
theta=linspace(0,90,nthe);
phi0=linspace(0,360,nphi);
phi=matrix(1,nphi);
E2sp=matrix(length(theta),length(phi0),length(f),6);
for(ii=1:length(phi0)) {
  phi(1,ii)=phi0(ii);
for(nf=1:length(f)) {
  E2sp(1:length(theta),ii,nf,1)=farfieldspherical( pinch(E2(1:na,1:nb,nf)), ux, uy, theta, phi0(ii));
}
}

#The mirror symmetry is used to calculate the far field field from dipoles located in region 1"  
for(nf=1:length(f)) {

E2sp_ms=matrix(length(theta),nphi);
i60=find(phi0,60);
i300=find(phi0,300);
for(ii=1:i60) {
  E2sp_ms(1:length(theta),i60-(ii-1))=E2sp(1:length(theta),ii,nf,1);
 }
for(ii=i60+1:length(phi0)) {
  E2sp_ms(1:length(theta),length(phi0)-(ii-i60))=E2sp(1:length(theta),ii,nf,1);
 }

#sum of the far fields from regions 1' and 1"
E2sp(1:length(theta),1:length(phi0),nf,1)=E2sp_ms(1:length(theta),1:length(phi0))+pinch(E2sp(1:length(theta),1:length(phi0),nf,1));
}

#60 rotational symmetry is used to get the farfields from dipoles located in regions 2 to 6 
for(nf=1:length(f)) {

for(jj=2:6) {
  i60=find(phi0,60);
  i300=find(phi0,300);
  for(ii=1:i300) {
    E2sp(1:nthe,(ii-1)+i60,nf,jj)=E2sp(1:nthe,ii,nf,jj-1);
     }
  for(ii=i300:length(phi0)-1) {
    E2sp(1:nthe,ii+1-i300,nf,jj)=E2sp(1:nthe,ii,nf,jj-1);
    }
}
}

#sum of the far fields from regions 1 to 6
for(nf=1:length(f)) {

for(ii=1:6) {
E2spA_tot(1:nthe,1:nphi,nf)=pinch(E2sp(1:nthe,1:nphi,nf,ii))+pinch(E2spA_tot(1:nthe,1:nphi,nf));
}
}

} # nn
} #if(ireg_1p)
##### calculation for regions from 2 to 6 -end- 

##### calculation for regions A,B,C and D -start-
E2spB_tot=matrix(nthe,nphi,length(f));
if(ireg_Ap) {
for(nn=1:ireg_Ap) {

load("OLED_hex_simple_"+"Ap"+num2str(nn)+"_x");
?"Ap, x, nn="+num2str(nn);
E2_far=getresult("far_field_change_index","E2_far");
save;
ux=E2_far.ux;
uy=E2_far.uy;
E2x=E2_far.E2_far1;


load("OLED_hex_simple_"+"Ap"+num2str(nn)+"_y");
?"Ap, y, nn="+num2str(nn);
E2_far=getresult("far_field_change_index","E2_far");
save;
E2y=E2_far.E2_far1;

load("OLED_hex_simple_"+"Ap"+num2str(nn)+"_z");
?"Ap, z, nn="+num2str(nn);
E2_far=getresult("far_field_change_index","E2_far");
save;
E2z=E2_far.E2_far1;

E2=E2x+E2y+E2z;

#transfomation of far field coordinate from (ux,uy) to (theta,phi) using farfieldspherical function
#E2sp=matrix(length(theta),length(phi0),6);
theta=linspace(0,90,nthe);
phi0=linspace(0,360,nphi);
phi=matrix(1,nphi);
E2sp=matrix(length(theta),length(phi0),length(f),6);

for(ii=1:length(phi0)) {
  phi(1,ii)=phi0(ii);
  E2sp(1:length(theta),ii,1:length(f),1)=farfieldspherical(E2, ux, uy, theta, phi0(ii));
}

#The mirror symmetry is used to calculate the far field field from dipoles located in region A"  
for(nf=1:length(f)) {

E2sp_ms=matrix(length(theta),nphi);
i60=find(phi0,60);
i300=find(phi0,300);
for(ii=1:i60) {
  E2sp_ms(1:length(theta),i60-(ii-1))=E2sp(1:length(theta),ii,nf,1);
 }
for(ii=i60+1:length(phi0)) {
  E2sp_ms(1:length(theta),length(phi0)-(ii-i60))=E2sp(1:length(theta),ii,nf,1);
 }

#sum of the far fields from regions A' and A"
E2sp(1:length(theta),1:length(phi0),nf,1)=E2sp_ms(1:length(theta),1:length(phi0))+pinch(E2sp(1:length(theta),1:length(phi0),nf,1));
}

#The 60 degree rotational symmetry is used to get the farfields from dipoles located in regions B, C and D 
for(nf=1:length(f)) {

for(jj=2:6) {
  i60=find(phi0,60);
  i300=find(phi0,300);
  for(ii=1:i300) {
    E2sp(1:nthe,(ii-1)+i60,nf,jj)=E2sp(1:nthe,ii,nf,jj-1);
   }
  for(ii=i300:length(phi0)-1) {
    E2sp(1:nthe,ii+1-i300,nf,jj)=E2sp(1:nthe,ii,nf,jj-1);
   }
}
}

#sum of the avaraged far fields between unit cell 1 and unit cell 2
for(nf=1:length(f)) {

E2spB_tot(1:nthe,1:nphi,nf)=pinch((E2sp(1:nthe,1:nphi,nf,1)+E2sp(1:nthe,1:nphi,nf,3)+E2sp(1:nthe,1:nphi,nf,4)+E2sp(1:nthe,1:nphi,nf,6)))/2
+pinch(E2spB_tot(1:nthe,1:nphi,nf));
}

} # nn
} #if(iregB)
##### calculation for regions A,B,C and D -end-

E2sp_tot=E2spA_tot+E2spB_tot; #far field from OLED with PC

#Radiated power calcualated from far field pattern
power_radiated=matrix(length(f));
dphi=(phi0(2)-phi0(1))*pi/180;
dtheta=(theta(2)-theta(1))*pi/180;
S=matrix(length(f));
for(nf=1:length(f)) {
  S(nf)=0;
  for(ip=1:length(phi0)) { 
    if(ip==1 or ip==length(phi0)) {DP=dphi/2;}
    else {DP=dphi;}
    for(it=1:length(theta)) {
      if(it==1 or it==length(theta)) {DT=dtheta/2;}
      else {DT=dtheta;}
      S(nf) =S(nf) +E2sp_tot(it,ip,nf)*sin(theta(it)*pi/180)*DP*DT;
     } #it
  } #ip
} # nf

for(nf=1:length(f)) {
  ?power_radiated(nf)=0.5*sqrt(eps0*(1.5^2)/mu0)*S(nf);
} # nf
##### calculation of far field rafiated from dipoles distributed randomly in the unit cell with PC -end-


##### calculation of far field pattern and radiated power from OLED with no PC -start-
pnts=ireg_1p+ireg_Ap;
dp0x=matrix(length(f));
Pw0=matrix(length(f));
E2_x0=matrix(na,nb,length(f));
E2_y0=matrix(na,nb,length(f));
E2_z0=matrix(na,nb,length(f));
E2_xyz0=matrix(na,nb,length(f));

load("OLED_hex_simple0_x");
Pwx0=matrix(length(f));
E2_far=getresult("far_field_change_index","E2_far");
T_far=getresult("far_field_change_index","T_far");
save;
E2_x0=E2_far.E2_far1; #farfield from x-directed dipole
Pwx0=T_far.T_far1*sourcepower(f); #radiated power from x-directed dipole
for(nf=1:length(f)) { #transpose of result from x-directed dipole to get result from y-directed dipole
E2_y0(1:na,1:nb,nf)=transpose(pinch(E2_x0(1:na,1:nb,nf)));
}

load("OLED_hex_simple0_z");
Pwz0=matrix(length(f));
E2_far=getresult("far_field_change_index","E2_far");
T_far=getresult("far_field_change_index","T_far");
save;
E2_z0=E2_far.E2_far1; #farfield by z-directed dipole
Pwz0=T_far.T_far1*sourcepower(f); #radiated power from z-directed dipole

E2_xyz0=E2_x0+E2_y0+E2_z0; #far field from OLED without PC
##### calculation of far field pattern and radiated power from OLED with no PC -end-

power_radiated0=matrix(length(f));
enhancement=matrix(length(f));
for(nf=1:length(f)) {
  ?power_radiated0(nf)=(Pwx0(nf)*2+Pwz0(nf))*(ireg_1p*2*6+ireg_Ap*2*2);
  ?enhancement(nf)=power_radiated(nf)/power_radiated0(nf); #extraction enhancement
} # nf

lmd=matrix(nf);
lmd(1:1:nf)=c/f(nf:-1:1);
plot(lmd*1e9,real(enhancement(nf:-1:1)),"Wavelength [nm]","Enhancement");

# save data to .mat file if plotting in Matlab
matlabsave("L2b",theta,phi,E2spA_tot,E2spB_tot,E2sp_tot);

# Interpolate from spherical to direction cosine coorinates
# and plot data for a specified wavelength

# Variables used for interpolation
lambda0 = 400e-9; # wavelength to plot
res=151;
ux=linspace(-1,1,res);
uy=ux;
E2_uxuy = matrix(res,res);
index=find(f,c/lambda0);

for (k=1:res) {
  for (j=1:res) {

    # transformation
    theta_tmp = real(acos( sqrt(1-ux(k)^2-uy(j)^2) ));
    phi_tmp   = real(atan( uy(j)/ux(k) ));

    if ( (ux(k)< 0)             ) { phi_tmp = phi_tmp+pi;    }
    if ( (ux(k)>=0) & (uy(j)<0) ) { phi_tmp = phi_tmp+2*pi;  }

    theta_index = find(theta*pi/180,theta_tmp);
    phi_index   = find(phi*pi/180,  phi_tmp);
    
    E2_uxuy(k,j) = E2sp_tot(theta_index,phi_index,index);    
  }
}
image(ux,uy,E2_uxuy(1:res,1:res),"","","lambda0 = "+num2str(lambda0),"polar");

if (save_data ==1) {
savedata("analysis_data");
}