clear; closeall;

load("usr_reverse_time_1.fsp");

if (havedata("R","Ex")==0) {
  selectpartial("import_source");
  set("enabled",0);
  selectpartial("beam_source");
  set("enabled",1);

  run;
}

E_R = getresult("R","E");
H_R = getresult("R","H");
E_T = getresult("T","E");
H_T = getresult("T","H");


# transform fields for reverse time simulation.  Define all quantities required 
# to define a custom field profile
# The conjugate(E) and -conjugate(H) operations are required to reverse the time.  
# The additional factor of -1(H) is required to correct for the direction of 
# propagation of this source (negative Z direction)
x=E_R.x; y=E_R.y; z=E_R.z;
Ex = conj(E_R.Ex);
Ey = conj(E_R.Ey);
Ez = conj(E_R.Ez);
Hx = -1*-conj(H_R.Hx);
Hy = -1*-conj(H_R.Hy);
Hz = -1*-conj(H_R.Hz);

include_H=1;
axis1=1;
axis2=2;
min1=min(x); max1=max(x);
min2=min(y); max2=max(y);
depth = z;
wavelength = c / E_R.f;  # source wavelength
index = 1;   # specify background refractive index
n1 = length(x);
n2 = length(y);
power = abs(transmission("R"));  # source power
sourcefilename = "R";

usr_create_fld;



x=E_T.x; y=E_T.y; z=E_T.z;
Ex = conj(E_T.Ex);
Ey = conj(E_T.Ey);
Ez = conj(E_T.Ez);
Hx = -conj(H_T.Hx);
Hy = -conj(H_T.Hy);
Hz = -conj(H_T.Hz);

include_H=1;
axis1=1;
axis2=2;
min1=min(x); max1=max(x);
min2=min(y); max2=max(y);
depth = z;
wavelength = c / E_T.f;  # source wavelength
index = 2;   # specify background refractive index
n1 = length(x);
n2 = length(y);
power = abs(transmission("T"));  # source power
sourcefilename = "T";

usr_create_fld;




# import into 2nd simulation;
switchtolayout;
save("usr_reverse_time_2");


if (havedata("R","Ex")==0) {
  selectpartial("beam_source");
  set("enabled",0);
  selectpartial("import_source");
  set("enabled",1);

  # manually import source profiles here
  asapimport("import_source_T","T.fld");
  asapimport("import_source_R","R.fld");
  # then run simulation
  run;
}