# This simple example illustrate how to simulating a polymer 
# that changes index as is exposed to a beam of light.  In this 
# example, the polymer's refractive index change is proportional
# to the power absorption.
#
# The initial index is 1.5+0.01i.  
# The maximum index is 1.6+0.01i.  
# We assume the imaginary part of the index is constant. 

closeall; clear;

# define initial refractive index profile.
res=50;
n_imag      = 0.01i;
n_real_start= 1.5;
n_real_max  = 1.6;
x_data = linspace(-0.5e-6,6e-6,res);
y_data = linspace(-0.5e-6,6e-6,res);
z_data = linspace(-8e-6,12e-6,res);
n=matrix(res,res,res) + n_real_start + n_imag;
curing_factor=30e-4;  # conversion factor from absorbed power to delta_n. Units of m^3/W/s

# define for loop
npts=5;
t_start=0;
t_end=1;
t=linspace(t_start,t_end,npts);
dt=t(2)-t(1);

for (i=1:npts) {

  # load refractive index profile
  switchtolayout;
  select("polymer");
  importnk2(n,x_data,y_data,z_data);    

  # run simulation
  run(3);

  # get simulation data
  runanalysis("Pabs");
  m    = "Pabs";
  x    = getdata(m,"x");
  y    = getdata(m,"y");
  z    = getdata(m,"z");
  f    = getdata(m,"f");
  E2   = pinch( getdata(m,"E2")  );
  index= pinch( getdata(m,"n")    );
  Pabs = pinch( getdata(m,"Pabs") );
  Pabs = Pabs*sourcepower(f);  # convert to Watts/um^2

  # plot results
  image(x*1e6,z*1e6,pinch(E2,2,find(y,0)),"x (um)","y (um)","|E|^2, t="+num2str(t(i)));
  setplot("colorbar min",0);  setplot("colorbar max",2);
  exportfigure("materials_polymer_field"+num2str(i));
  image(x*1e6,z*1e6,pinch(Pabs,2,find(y,0)),"x (um)","y (um)","Absorption, t="+num2str(t(i)));
  setplot("colorbar min",0);  setplot("colorbar max",1500);
  exportfigure("materials_polymer_absorption"+num2str(i));
  image(x*1e6,z*1e6,real(pinch(index,2,find(y,0))),"x (um)","y (um)","Re(index), t="+num2str(t(i))); 
  setplot("colorbar min",1.4);  setplot("colorbar max",n_real_max);
  exportfigure("materials_polymer_index"+num2str(i));
  ?"T="+num2str(transmission("T"))+",  R="+num2str(-transmission("R"));

  # calculate new refractive index profile
  # for next pass throught the loop
  Pabs=interp(Pabs,x,y,z,x_data,y_data,z_data);                  # interpolate from simulation grid to data grid
  n     = real(n);                                      # take real part
  n     = n + Pabs*dt*curing_factor;                                 # increase refractive index by some factor of the absorption
  n     = (n<n_real_max)*n + (n>n_real_max)*n_real_max; # don't allow any values over n_real_max
  n     = n + n_imag;                                   # add imaginary part back in

}